From 7030361aab7e9d7927d98a5c68211c05531b256d Mon Sep 17 00:00:00 2001 From: Berk Hess Date: Tue, 10 Dec 2019 10:56:19 +0100 Subject: [PATCH] Use ListOfLists for atom to constraint mapping Replace t_blocka by ListOfList for all atom to constraint map handling. Change-Id: If33e7068949c28453a186c4a83842bab3dc6fb2a --- src/gromacs/domdec/domdec_constraints.cpp | 72 ++++---- src/gromacs/gmxpreprocess/grompp.cpp | 21 +-- src/gromacs/mdlib/constr.cpp | 65 ++++--- src/gromacs/mdlib/constr.h | 44 ++--- src/gromacs/mdlib/constraintrange.cpp | 33 ++-- src/gromacs/mdlib/lincs.cpp | 164 ++++++++---------- src/gromacs/mdlib/lincs.h | 20 +-- src/gromacs/mdlib/tests/constrtestrunners.cpp | 9 +- src/gromacs/mdlib/updategroups.cpp | 52 +++--- 9 files changed, 218 insertions(+), 262 deletions(-) diff --git a/src/gromacs/domdec/domdec_constraints.cpp b/src/gromacs/domdec/domdec_constraints.cpp index 0a62ae980e..99b79d828d 100644 --- a/src/gromacs/domdec/domdec_constraints.cpp +++ b/src/gromacs/domdec/domdec_constraints.cpp @@ -67,11 +67,13 @@ #include "gromacs/utility/exceptions.h" #include "gromacs/utility/fatalerror.h" #include "gromacs/utility/gmxassert.h" -#include "gromacs/utility/smalloc.h" +#include "gromacs/utility/listoflists.h" #include "domdec_internal.h" #include "domdec_specatomcomm.h" +using gmx::ListOfLists; + /*! \brief Struct used during constraint setup with domain decomposition */ struct gmx_domdec_constraints_t { @@ -141,7 +143,7 @@ static void walk_out(int con, int nrec, gmx::ArrayRef ia1, gmx::ArrayRef ia2, - const t_blocka* at2con, + const ListOfLists& at2con, const gmx_ga2la_t& ga2la, gmx_bool bHomeConnect, gmx_domdec_constraints_t* dc, @@ -149,9 +151,6 @@ static void walk_out(int con, t_ilist* il_local, std::vector* ireq) { - int a1_gl, a2_gl, i, coni, b; - const t_iatom* iap; - if (!dc->gc_req[con_offset + con]) { /* Add this non-home constraint to the list */ @@ -163,10 +162,10 @@ static void walk_out(int con, il_local->nalloc = over_alloc_dd(il_local->nr + 3); srenew(il_local->iatoms, il_local->nalloc); } - iap = constr_iatomptr(ia1, ia2, con); + const int* iap = constr_iatomptr(ia1, ia2, con); il_local->iatoms[il_local->nr++] = iap[0]; - a1_gl = offset + iap[1]; - a2_gl = offset + iap[2]; + const int a1_gl = offset + iap[1]; + const int a2_gl = offset + iap[2]; /* The following indexing code can probably be optizimed */ if (const int* a_loc = ga2la.findHome(a1_gl)) { @@ -200,13 +199,14 @@ static void walk_out(int con, if (nrec > 0) { - for (i = at2con->index[a]; i < at2con->index[a + 1]; i++) + /* Loop over the constraint connected to atom a */ + for (const int coni : at2con[a]) { - coni = at2con->a[i]; if (coni != con) { /* Walk further */ - iap = constr_iatomptr(ia1, ia2, coni); + const int* iap = constr_iatomptr(ia1, ia2, coni); + int b; if (a == iap[1]) { b = iap[2]; @@ -306,17 +306,14 @@ static void atoms_to_settles(gmx_domdec_t* dd, } /*! \brief Looks up constraint for the local atoms */ -static void atoms_to_constraints(gmx_domdec_t* dd, - const gmx_mtop_t* mtop, - const int* cginfo, - gmx::ArrayRef at2con_mt, - int nrec, - t_ilist* ilc_local, - std::vector* ireq) +static void atoms_to_constraints(gmx_domdec_t* dd, + const gmx_mtop_t* mtop, + const int* cginfo, + gmx::ArrayRef> at2con_mt, + int nrec, + t_ilist* ilc_local, + std::vector* ireq) { - const t_blocka* at2con; - int b_lo, offset, b_mol, i, con, con_offset; - gmx_domdec_constraints_t* dc = dd->constraints; gmx_domdec_specat_comm_t* dcc = dd->constraint_comm; @@ -344,15 +341,16 @@ static void atoms_to_constraints(gmx_domdec_t* dd, * This is only required for the global index to make sure * that we use each constraint only once. */ - con_offset = dc->molb_con_offset[mb] + molnr * dc->molb_ncon_mol[mb]; + const int con_offset = dc->molb_con_offset[mb] + molnr * dc->molb_ncon_mol[mb]; /* The global atom number offset for this molecule */ - offset = a_gl - a_mol; - at2con = &at2con_mt[molb.type]; - for (i = at2con->index[a_mol]; i < at2con->index[a_mol + 1]; i++) + const int offset = a_gl - a_mol; + /* Loop over the constraints connected to atom a_mol in the molecule */ + const auto& at2con = at2con_mt[molb.type]; + for (const int con : at2con[a_mol]) { - con = at2con->a[i]; const int* iap = constr_iatomptr(ia1, ia2, con); + int b_mol; if (a_mol == iap[1]) { b_mol = iap[2]; @@ -373,7 +371,7 @@ static void atoms_to_constraints(gmx_domdec_t* dd, ilc_local->nalloc = over_alloc_dd(ilc_local->nr + 3); srenew(ilc_local->iatoms, ilc_local->nalloc); } - b_lo = *a_loc; + 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); @@ -416,13 +414,11 @@ int dd_make_local_constraints(gmx_domdec_t* dd, int nrec, t_ilist* il_local) { - gmx_domdec_constraints_t* dc; - t_ilist * ilc_local, *ils_local; - std::vector* ireq; - gmx::ArrayRef at2con_mt; - gmx::HashedMap* ga2la_specat; - int at_end, i, j; - t_iatom* iap; + gmx_domdec_constraints_t* dc; + t_ilist * 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 @@ -446,6 +442,8 @@ int dd_make_local_constraints(gmx_domdec_t* dd, dc->ncon = 0; ilc_local->nr = 0; + gmx::ArrayRef> at2con_mt; + std::vector* ireq = nullptr; if (dd->constraint_comm) { // TODO Perhaps gmx_domdec_constraints_t should keep a valid constr? @@ -454,12 +452,6 @@ int dd_make_local_constraints(gmx_domdec_t* dd, ireq = &dc->requestedGlobalAtomIndices[0]; ireq->clear(); } - else - { - // Currently unreachable - at2con_mt = {}; - ireq = nullptr; - } gmx::ArrayRef> at2settle_mt; /* When settle works inside charge groups, we assigned them already */ diff --git a/src/gromacs/gmxpreprocess/grompp.cpp b/src/gromacs/gmxpreprocess/grompp.cpp index 735b5c37d6..21f967e02a 100644 --- a/src/gromacs/gmxpreprocess/grompp.cpp +++ b/src/gromacs/gmxpreprocess/grompp.cpp @@ -102,6 +102,7 @@ #include "gromacs/utility/futil.h" #include "gromacs/utility/gmxassert.h" #include "gromacs/utility/keyvaluetreebuilder.h" +#include "gromacs/utility/listoflists.h" #include "gromacs/utility/mdmodulenotification.h" #include "gromacs/utility/smalloc.h" #include "gromacs/utility/snprintf.h" @@ -1426,7 +1427,7 @@ static bool haveDecoupledModeInMol(const gmx_moltype_t& molt, const t_atom* atom = molt.atoms.atom; - t_blocka atomToConstraints = + const auto atomToConstraints = gmx::make_at2con(molt, iparams, gmx::FlexibleConstraintTreatment::Exclude); bool haveDecoupledMode = false; @@ -1456,23 +1457,21 @@ static bool haveDecoupledModeInMol(const gmx_moltype_t& molt, int a1 = il.iatoms[1 + i + 1]; int a2 = il.iatoms[1 + i + 2]; if ((atom[a0].m > atom[a2].m * massFactorThreshold || atom[a2].m > atom[a0].m * massFactorThreshold) - && atomToConstraints.index[a0 + 1] - atomToConstraints.index[a0] == 1 - && atomToConstraints.index[a2 + 1] - atomToConstraints.index[a2] == 1 - && atomToConstraints.index[a1 + 1] - atomToConstraints.index[a1] >= 3) + && atomToConstraints[a0].ssize() == 1 && atomToConstraints[a2].ssize() == 1 + && atomToConstraints[a1].ssize() >= 3) { - int constraint0 = atomToConstraints.a[atomToConstraints.index[a0]]; - int constraint2 = atomToConstraints.a[atomToConstraints.index[a2]]; + int constraint0 = atomToConstraints[a0][0]; + int constraint2 = atomToConstraints[a2][0]; bool foundAtom0 = false; bool foundAtom2 = false; - for (int conIndex = atomToConstraints.index[a1]; - conIndex < atomToConstraints.index[a1 + 1]; conIndex++) + for (const int constraint : atomToConstraints[a1]) { - if (atomToConstraints.a[conIndex] == constraint0) + if (constraint == constraint0) { foundAtom0 = true; } - if (atomToConstraints.a[conIndex] == constraint2) + if (constraint == constraint2) { foundAtom2 = true; } @@ -1486,8 +1485,6 @@ static bool haveDecoupledModeInMol(const gmx_moltype_t& molt, } } - done_blocka(&atomToConstraints); - return haveDecoupledMode; } diff --git a/src/gromacs/mdlib/constr.cpp b/src/gromacs/mdlib/constr.cpp index 984c901e62..1e418cc78e 100644 --- a/src/gromacs/mdlib/constr.cpp +++ b/src/gromacs/mdlib/constr.cpp @@ -72,7 +72,6 @@ #include "gromacs/pbcutil/pbc.h" #include "gromacs/pulling/pull.h" #include "gromacs/timing/wallcycle.h" -#include "gromacs/topology/block.h" #include "gromacs/topology/ifunc.h" #include "gromacs/topology/mtop_lookup.h" #include "gromacs/topology/mtop_util.h" @@ -80,6 +79,7 @@ #include "gromacs/utility/exceptions.h" #include "gromacs/utility/fatalerror.h" #include "gromacs/utility/gmxassert.h" +#include "gromacs/utility/listoflists.h" #include "gromacs/utility/pleasecite.h" #include "gromacs/utility/smalloc.h" #include "gromacs/utility/txtdump.h" @@ -132,7 +132,7 @@ public: //! The number of flexible constraints. int nflexcon = 0; //! A list of atoms to constraints for each moleculetype. - std::vector at2con_mt; + std::vector> at2con_mt; //! A list of atoms to settles for each moleculetype std::vector> at2settle_mt; //! LINCS data. @@ -743,14 +743,17 @@ FlexibleConstraintTreatment flexibleConstraintTreatment(bool haveDynamicsIntegra * \param[in] numAtoms The number of atoms to construct the list for * \param[in] ilists The interaction lists, size F_NRE * \param[in] iparams Interaction parameters, can be null when - * flexibleConstraintTreatment=Include \param[in] flexibleConstraintTreatment The flexible - * constraint treatment, see enum above \returns a block struct with all constraints for each atom + * \p flexibleConstraintTreatment==Include + * \param[in] flexibleConstraintTreatment The flexible constraint treatment, + * see enum above + * + * \returns a block struct with all constraints for each atom */ template -static t_blocka makeAtomsToConstraintsList(int numAtoms, - const T* ilists, - const t_iparams* iparams, - FlexibleConstraintTreatment flexibleConstraintTreatment) +static ListOfLists makeAtomsToConstraintsList(int numAtoms, + const T* ilists, + const t_iparams* iparams, + FlexibleConstraintTreatment flexibleConstraintTreatment) { GMX_ASSERT(flexibleConstraintTreatment == FlexibleConstraintTreatment::Include || iparams != nullptr, "With flexible constraint detection we need valid iparams"); @@ -775,19 +778,13 @@ static t_blocka makeAtomsToConstraintsList(int numAtoms, } } - t_blocka at2con; - at2con.nr = numAtoms; - at2con.nalloc_index = at2con.nr + 1; - snew(at2con.index, at2con.nalloc_index); - at2con.index[0] = 0; + std::vector listRanges(numAtoms + 1); for (int a = 0; a < numAtoms; a++) { - at2con.index[a + 1] = at2con.index[a] + count[a]; - count[a] = 0; + listRanges[a + 1] = listRanges[a] + count[a]; + count[a] = 0; } - at2con.nra = at2con.index[at2con.nr]; - at2con.nalloc_a = at2con.nra; - snew(at2con.a, at2con.nalloc_a); + std::vector elements(listRanges[numAtoms]); /* The F_CONSTRNC constraints have constraint numbers * that continue after the last F_CONSTR constraint. @@ -804,28 +801,28 @@ static t_blocka makeAtomsToConstraintsList(int numAtoms, { for (int j = 1; j < 3; j++) { - int a = ilist.iatoms[i + j]; - at2con.a[at2con.index[a] + count[a]++] = numConstraints; + const int a = ilist.iatoms[i + j]; + elements[listRanges[a] + count[a]++] = numConstraints; } } numConstraints++; } } - return at2con; + return ListOfLists(std::move(listRanges), std::move(elements)); } -t_blocka make_at2con(int numAtoms, - const t_ilist* ilist, - const t_iparams* iparams, - FlexibleConstraintTreatment flexibleConstraintTreatment) +ListOfLists make_at2con(int numAtoms, + const t_ilist* ilist, + const t_iparams* iparams, + FlexibleConstraintTreatment flexibleConstraintTreatment) { return makeAtomsToConstraintsList(numAtoms, ilist, iparams, flexibleConstraintTreatment); } -t_blocka make_at2con(const gmx_moltype_t& moltype, - gmx::ArrayRef iparams, - FlexibleConstraintTreatment flexibleConstraintTreatment) +ListOfLists make_at2con(const gmx_moltype_t& moltype, + gmx::ArrayRef iparams, + FlexibleConstraintTreatment flexibleConstraintTreatment) { return makeAtomsToConstraintsList(moltype.atoms.nr, moltype.ilist.data(), iparams.data(), flexibleConstraintTreatment); @@ -924,10 +921,10 @@ void Constraints::setConstraints(const gmx_localtop_t& top, const t_mdatoms& md) * indices to constraint indices. * * Note that flexible constraints are only enabled with a dynamical integrator. */ -static std::vector makeAtomToConstraintMappings(const gmx_mtop_t& mtop, - FlexibleConstraintTreatment flexibleConstraintTreatment) +static std::vector> makeAtomToConstraintMappings(const gmx_mtop_t& mtop, + FlexibleConstraintTreatment flexibleConstraintTreatment) { - std::vector mapping; + std::vector> mapping; mapping.reserve(mtop.moltype.size()); for (const gmx_moltype_t& moltype : mtop.moltype) { @@ -1094,10 +1091,6 @@ Constraints::Impl::Impl(const gmx_mtop_t& mtop_p, Constraints::Impl::~Impl() { - for (auto blocka : at2con_mt) - { - done_blocka(&blocka); - } if (bSettleErrorHasOccurred != nullptr) { sfree(bSettleErrorHasOccurred); @@ -1118,7 +1111,7 @@ void Constraints::saveEdsamPointer(gmx_edsam* ed) impl_->ed = ed; } -ArrayRef Constraints::atom2constraints_moltype() const +ArrayRef> Constraints::atom2constraints_moltype() const { return impl_->at2con_mt; } diff --git a/src/gromacs/mdlib/constr.h b/src/gromacs/mdlib/constr.h index 2ddc2704d1..a6d159f92a 100644 --- a/src/gromacs/mdlib/constr.h +++ b/src/gromacs/mdlib/constr.h @@ -63,7 +63,6 @@ struct gmx_mtop_t; struct gmx_multisim_t; struct gmx_wallcycle; struct pull_t; -struct t_blocka; struct t_commrec; struct t_ilist; struct t_inputrec; @@ -76,6 +75,8 @@ namespace gmx { template class ArrayRefWithPadding; +template +class ListOfLists; //! Describes supported flavours of constrained updates. enum class ConstraintVariable : int @@ -183,7 +184,7 @@ public: //! Links the essentialdynamics and constraint code. void saveEdsamPointer(gmx_edsam* ed); //! Getter for use by domain decomposition. - ArrayRef atom2constraints_moltype() const; + ArrayRef> atom2constraints_moltype() const; //! Getter for use by domain decomposition. ArrayRef> atom2settle_moltype() const; @@ -234,40 +235,43 @@ enum class FlexibleConstraintTreatment /*! \brief Returns the flexible constraint treatment depending on whether the integrator is dynamic */ FlexibleConstraintTreatment flexibleConstraintTreatment(bool haveDynamicsIntegrator); -/*! \brief Returns a block struct to go from atoms to constraints +/*! \brief Returns a ListOfLists object to go from atoms to constraints * - * The block struct will contain constraint indices with lower indices + * The object will contain constraint indices with lower indices * directly matching the order in F_CONSTR and higher indices matching * the order in F_CONSTRNC offset by the number of constraints in F_CONSTR. * * \param[in] moltype The molecule data * \param[in] iparams Interaction parameters, can be null when - * flexibleConstraintTreatment=Include \param[in] flexibleConstraintTreatment The flexible - * constraint treatment, see enum above \returns a block struct with all constraints for each atom + * \p flexibleConstraintTreatment==Include + * \param[in] flexibleConstraintTreatment The flexible constraint treatment, + * see enum above + * + * \returns a ListOfLists object with all constraints for each atom */ -t_blocka make_at2con(const gmx_moltype_t& moltype, - gmx::ArrayRef iparams, - FlexibleConstraintTreatment flexibleConstraintTreatment); +ListOfLists make_at2con(const gmx_moltype_t& moltype, + gmx::ArrayRef iparams, + FlexibleConstraintTreatment flexibleConstraintTreatment); -/*! \brief Returns a block struct to go from atoms to constraints +/*! \brief Returns a ListOfLists object to go from atoms to constraints * - * The block struct will contain constraint indices with lower indices + * The object will contain constraint indices with lower indices * directly matching the order in F_CONSTR and higher indices matching * the order in F_CONSTRNC offset by the number of constraints in F_CONSTR. * * \param[in] numAtoms The number of atoms to construct the list for * \param[in] ilist Interaction list, size F_NRE * \param[in] iparams Interaction parameters, can be null when - * flexibleConstraintTreatment=Include \param[in] flexibleConstraintTreatment The flexible - * constraint treatment, see enum above \returns a block struct with all constraints for each atom + * \p flexibleConstraintTreatment==Include + * \param[in] flexibleConstraintTreatment The flexible constraint treatment, + * see enum above + * + * \returns a ListOfLists object with all constraints for each atom */ -t_blocka make_at2con(int numAtoms, - const t_ilist* ilist, - const t_iparams* iparams, - FlexibleConstraintTreatment flexibleConstraintTreatment); - -/*! \brief Returns an array of atom to constraints lists for the moltypes */ -const t_blocka* atom2constraints_moltype(const Constraints* constr); +ListOfLists make_at2con(int numAtoms, + const t_ilist* ilist, + const t_iparams* 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); diff --git a/src/gromacs/mdlib/constraintrange.cpp b/src/gromacs/mdlib/constraintrange.cpp index 7e5b320e4d..10f00fa29c 100644 --- a/src/gromacs/mdlib/constraintrange.cpp +++ b/src/gromacs/mdlib/constraintrange.cpp @@ -45,32 +45,30 @@ #include "gromacs/mdlib/constr.h" #include "gromacs/mdtypes/inputrec.h" -#include "gromacs/topology/block.h" #include "gromacs/topology/mtop_util.h" #include "gromacs/utility/basedefinitions.h" #include "gromacs/utility/fatalerror.h" +#include "gromacs/utility/listoflists.h" #include "gromacs/utility/logger.h" #include "gromacs/utility/real.h" -#include "gromacs/utility/smalloc.h" namespace gmx { //! Recursing function to help find all adjacent constraints. -static void constr_recur(const t_blocka* at2con, +static void constr_recur(const ListOfLists& at2con, const InteractionLists& ilist, gmx::ArrayRef iparams, gmx_bool bTopB, int at, int depth, int nc, - int* path, + ArrayRef path, real r0, real r1, real* r2max, int* count) { - int c, con, a1; gmx_bool bUse; real len, rn0, rn1; @@ -80,12 +78,11 @@ static void constr_recur(const t_blocka* at2con, gmx::ArrayRef ia2 = ilist[F_CONSTRNC].iatoms; /* Loop over all constraints connected to this atom */ - for (c = at2con->index[at]; c < at2con->index[at + 1]; c++) + for (const int con : at2con[at]) { - con = at2con->a[c]; /* Do not walk over already used constraints */ bUse = TRUE; - for (a1 = 0; a1 < depth; a1++) + for (int a1 = 0; a1 < depth; a1++) { if (con == path[a1]) { @@ -124,7 +121,7 @@ static void constr_recur(const t_blocka* at2con, fprintf(debug, "Found longer constraint distance: r0 %5.3f r1 %5.3f rmax %5.3f\n", rn0, rn1, sqrt(*r2max)); - for (a1 = 0; a1 < depth; a1++) + for (int a1 = 0; a1 < depth; a1++) { fprintf(debug, " %d %5.3f", path[a1], iparams[constr_iatomptr(ia1, ia2, con)[0]].constr.dA); @@ -138,6 +135,7 @@ static void constr_recur(const t_blocka* at2con, */ if (depth + 1 < nc && *count < 1000 * nc) { + int a1; if (ia[1] == at) { a1 = ia[2]; @@ -160,10 +158,9 @@ static real constr_r_max_moltype(const gmx_moltype_t* molt, gmx::ArrayRef iparams, const t_inputrec* ir) { - int natoms, *path, at, count; + int natoms, at, count; - t_blocka at2con; - real r0, r1, r2maxA, r2maxB, rmax, lam0, lam1; + real r0, r1, r2maxA, r2maxB, rmax, lam0, lam1; if (molt->ilist[F_CONSTR].size() == 0 && molt->ilist[F_CONSTRNC].size() == 0) { @@ -172,8 +169,9 @@ static real constr_r_max_moltype(const gmx_moltype_t* molt, natoms = molt->atoms.nr; - at2con = make_at2con(*molt, iparams, flexibleConstraintTreatment(EI_DYNAMICS(ir->eI))); - snew(path, 1 + ir->nProjOrder); + const ListOfLists at2con = + make_at2con(*molt, iparams, flexibleConstraintTreatment(EI_DYNAMICS(ir->eI))); + std::vector path(1 + ir->nProjOrder); for (at = 0; at < 1 + ir->nProjOrder; at++) { path[at] = -1; @@ -186,7 +184,7 @@ static real constr_r_max_moltype(const gmx_moltype_t* molt, r1 = 0; count = 0; - constr_recur(&at2con, molt->ilist, iparams, FALSE, at, 0, 1 + ir->nProjOrder, path, r0, r1, + constr_recur(at2con, molt->ilist, iparams, FALSE, at, 0, 1 + ir->nProjOrder, path, r0, r1, &r2maxA, &count); } if (ir->efep == efepNO) @@ -201,7 +199,7 @@ static real constr_r_max_moltype(const gmx_moltype_t* molt, r0 = 0; r1 = 0; count = 0; - constr_recur(&at2con, molt->ilist, iparams, TRUE, at, 0, 1 + ir->nProjOrder, path, r0, + constr_recur(at2con, molt->ilist, iparams, TRUE, at, 0, 1 + ir->nProjOrder, path, r0, r1, &r2maxB, &count); } lam0 = ir->fepvals->init_lambda; @@ -217,9 +215,6 @@ static real constr_r_max_moltype(const gmx_moltype_t* molt, } } - done_blocka(&at2con); - sfree(path); - return rmax; } diff --git a/src/gromacs/mdlib/lincs.cpp b/src/gromacs/mdlib/lincs.cpp index 81fce32355..9e4b18f304 100644 --- a/src/gromacs/mdlib/lincs.cpp +++ b/src/gromacs/mdlib/lincs.cpp @@ -73,7 +73,6 @@ #include "gromacs/simd/simd.h" #include "gromacs/simd/simd_math.h" #include "gromacs/simd/vector_operations.h" -#include "gromacs/topology/block.h" #include "gromacs/topology/mtop_util.h" #include "gromacs/utility/alignedallocator.h" #include "gromacs/utility/arrayref.h" @@ -83,6 +82,7 @@ #include "gromacs/utility/exceptions.h" #include "gromacs/utility/fatalerror.h" #include "gromacs/utility/gmxomp.h" +#include "gromacs/utility/listoflists.h" #include "gromacs/utility/pleasecite.h" using namespace gmx; // TODO: Remove when this file is moved into gmx namespace @@ -1342,33 +1342,29 @@ static void set_lincs_matrix(Lincs* li, real* invmass, real lambda) } //! Finds all triangles of atoms that share constraints to a central atom. -static int count_triangle_constraints(const InteractionLists& ilist, const t_blocka& at2con) +static int count_triangle_constraints(const InteractionLists& ilist, const ListOfLists& at2con) { - int ncon1, ncon_tot; - int c0, n1, c1, ac1, n2, c2; - int ncon_triangle; - - ncon1 = ilist[F_CONSTR].size() / 3; - ncon_tot = ncon1 + ilist[F_CONSTRNC].size() / 3; + const int ncon1 = ilist[F_CONSTR].size() / 3; + const int ncon_tot = ncon1 + ilist[F_CONSTRNC].size() / 3; gmx::ArrayRef ia1 = ilist[F_CONSTR].iatoms; gmx::ArrayRef ia2 = ilist[F_CONSTRNC].iatoms; - ncon_triangle = 0; - for (c0 = 0; c0 < ncon_tot; c0++) + int ncon_triangle = 0; + for (int c0 = 0; c0 < ncon_tot; c0++) { bool bTriangle = FALSE; const int* iap = constr_iatomptr(ia1, ia2, c0); const int a00 = iap[1]; const int a01 = iap[2]; - for (n1 = at2con.index[a01]; n1 < at2con.index[a01 + 1]; n1++) + for (const int c1 : at2con[a01]) { - c1 = at2con.a[n1]; if (c1 != c0) { const int* iap = constr_iatomptr(ia1, ia2, c1); const int a10 = iap[1]; const int a11 = iap[2]; + int ac1; if (a10 == a01) { ac1 = a11; @@ -1377,9 +1373,8 @@ static int count_triangle_constraints(const InteractionLists& ilist, const t_blo { ac1 = a10; } - for (n2 = at2con.index[ac1]; n2 < at2con.index[ac1 + 1]; n2++) + for (const int c2 : at2con[ac1]) { - c2 = at2con.a[n2]; if (c2 != c0 && c2 != c1) { const int* iap = constr_iatomptr(ia1, ia2, c2); @@ -1403,40 +1398,36 @@ static int count_triangle_constraints(const InteractionLists& ilist, const t_blo } //! Finds sequences of sequential constraints. -static bool more_than_two_sequential_constraints(const InteractionLists& ilist, const t_blocka& at2con) +static bool more_than_two_sequential_constraints(const InteractionLists& ilist, const ListOfLists& at2con) { - int ncon1, ncon_tot, c; - bool bMoreThanTwoSequentialConstraints; - - ncon1 = ilist[F_CONSTR].size() / 3; - ncon_tot = ncon1 + ilist[F_CONSTRNC].size() / 3; + const int ncon1 = ilist[F_CONSTR].size() / 3; + const int ncon_tot = ncon1 + ilist[F_CONSTRNC].size() / 3; gmx::ArrayRef ia1 = ilist[F_CONSTR].iatoms; gmx::ArrayRef ia2 = ilist[F_CONSTRNC].iatoms; - bMoreThanTwoSequentialConstraints = FALSE; - for (c = 0; c < ncon_tot && !bMoreThanTwoSequentialConstraints; c++) + for (int c = 0; c < ncon_tot; c++) { const int* iap = constr_iatomptr(ia1, ia2, c); const int a1 = iap[1]; const int a2 = iap[2]; /* Check if this constraint has constraints connected at both atoms */ - if (at2con.index[a1 + 1] - at2con.index[a1] > 1 && at2con.index[a2 + 1] - at2con.index[a2] > 1) + if (at2con[a1].ssize() > 1 && at2con[a2].ssize() > 1) { - bMoreThanTwoSequentialConstraints = TRUE; + return true; } } - return bMoreThanTwoSequentialConstraints; + return false; } -Lincs* init_lincs(FILE* fplog, - const gmx_mtop_t& mtop, - int nflexcon_global, - ArrayRef at2con, - bool bPLINCS, - int nIter, - int nProjOrder) +Lincs* init_lincs(FILE* fplog, + const gmx_mtop_t& mtop, + int nflexcon_global, + ArrayRef> atomToConstraintsPerMolType, + bool bPLINCS, + int nIter, + int nProjOrder) { // TODO this should become a unique_ptr Lincs* li; @@ -1458,9 +1449,10 @@ Lincs* init_lincs(FILE* fplog, li->max_connect = 0; for (size_t mt = 0; mt < mtop.moltype.size(); mt++) { + const auto& at2con = atomToConstraintsPerMolType[mt]; for (int a = 0; a < mtop.moltype[mt].atoms.nr; a++) { - li->max_connect = std::max(li->max_connect, at2con[mt].index[a + 1] - at2con[mt].index[a]); + li->max_connect = std::max(li->max_connect, int(at2con[a].ssize())); } } @@ -1468,11 +1460,12 @@ Lincs* init_lincs(FILE* fplog, bMoreThanTwoSeq = FALSE; for (const gmx_molblock_t& molb : mtop.molblock) { - const gmx_moltype_t& molt = mtop.moltype[molb.type]; + const gmx_moltype_t& molt = mtop.moltype[molb.type]; + const auto& at2con = atomToConstraintsPerMolType[molb.type]; - li->ncg_triangle += molb.nmol * count_triangle_constraints(molt.ilist, at2con[molb.type]); + li->ncg_triangle += molb.nmol * count_triangle_constraints(molt.ilist, at2con); - if (!bMoreThanTwoSeq && more_than_two_sequential_constraints(molt.ilist, at2con[molb.type])) + if (!bMoreThanTwoSeq && more_than_two_sequential_constraints(molt.ilist, at2con)) { bMoreThanTwoSeq = TRUE; } @@ -1645,7 +1638,13 @@ static void lincs_thread_setup(Lincs* li, int natoms) } //! Assign a constraint. -static void assign_constraint(Lincs* li, int constraint_index, int a1, int a2, real lenA, real lenB, const t_blocka* at2con) +static void assign_constraint(Lincs* li, + int constraint_index, + int a1, + int a2, + real lenA, + real lenB, + const ListOfLists& at2con) { int con; @@ -1664,8 +1663,7 @@ static void assign_constraint(Lincs* li, int constraint_index, int a1, int a2, r /* Make space in the constraint connection matrix for constraints * connected to both end of the current constraint. */ - li->ncc += at2con->index[a1 + 1] - at2con->index[a1] - 1 + at2con->index[a2 + 1] - - at2con->index[a2] - 1; + li->ncc += at2con[a1].ssize() - 1 + at2con[a2].ssize() - 1; li->blnr[con + 1] = li->ncc; @@ -1675,13 +1673,13 @@ static void assign_constraint(Lincs* li, int constraint_index, int a1, int a2, r /*! \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 t_blocka* at2con) +static void check_assign_connected(Lincs* li, + const t_iatom* iatom, + const t_idef& 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 @@ -1690,28 +1688,18 @@ static void check_assign_connected(Lincs* li, * connected constraints. We need to assign those * to the same task. */ - int end; - - for (end = 0; end < 2; end++) + for (int end = 0; end < 2; end++) { - int a, k; + const int a = (end == 0 ? a1 : a2); - a = (end == 0 ? a1 : a2); - - for (k = at2con->index[a]; k < at2con->index[a + 1]; k++) + for (const int cc : at2con[a]) { - int cc; - - cc = at2con->a[k]; /* Check if constraint cc has not yet been assigned */ if (li->con_index[cc] == -1) { - int type; - real lenA, lenB; - - type = iatom[cc * 3]; - lenA = idef.iparams[type].constr.dA; - lenB = idef.iparams[type].constr.dB; + const int type = iatom[cc * 3]; + const real lenA = idef.iparams[type].constr.dA; + const real lenB = idef.iparams[type].constr.dB; if (bDynamics || lenA != 0 || lenB != 0) { @@ -1725,24 +1713,21 @@ 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 t_blocka* at2con) +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) { - int nca, cc[32], ca[32], k; + int nca, cc[32], ca[32]; int c_triangle[2] = { -1, -1 }; nca = 0; - for (k = at2con->index[a1]; k < at2con->index[a1 + 1]; k++) + for (const int c : at2con[a1]) { - int c; - - c = at2con->a[k]; if (c != constraint_index) { int aa1, aa2; @@ -1764,11 +1749,8 @@ static void check_assign_triangle(Lincs* li, } } - for (k = at2con->index[a2]; k < at2con->index[a2 + 1]; k++) + for (const int c : at2con[a2]) { - int c; - - c = at2con->a[k]; if (c != constraint_index) { int aa1, aa2, i; @@ -1827,7 +1809,7 @@ static void check_assign_triangle(Lincs* li, } //! Sets matrix indices. -static void set_matrix_indices(Lincs* li, const Task& li_task, const t_blocka* at2con, bool bSortMatrix) +static void set_matrix_indices(Lincs* li, const Task& li_task, const ListOfLists& at2con, bool bSortMatrix) { for (int b = li_task.b0; b < li_task.b1; b++) { @@ -1835,17 +1817,17 @@ static void set_matrix_indices(Lincs* li, const Task& li_task, const t_blocka* a const int a2 = li->atoms[b].index2; int i = li->blnr[b]; - for (int k = at2con->index[a1]; k < at2con->index[a1 + 1]; k++) + for (const int constraint : at2con[a1]) { - int concon = li->con_index[at2con->a[k]]; + const int concon = li->con_index[constraint]; if (concon != b) { li->blbnb[i++] = concon; } } - for (int k = at2con->index[a2]; k < at2con->index[a2 + 1]; k++) + for (const int constraint : at2con[a2]) { - int concon = li->con_index[at2con->a[k]]; + const int concon = li->con_index[constraint]; if (concon != b) { li->blbnb[i++] = concon; @@ -1863,7 +1845,6 @@ static void set_matrix_indices(Lincs* li, const Task& li_task, const t_blocka* a void set_lincs(const t_idef& idef, const t_mdatoms& md, bool bDynamics, const t_commrec* cr, Lincs* li) { int natoms; - t_blocka at2con; t_iatom* iatom; li->nc_real = 0; @@ -1915,7 +1896,8 @@ void set_lincs(const t_idef& idef, const t_mdatoms& md, bool bDynamics, const t_ natoms = md.homenr; } - at2con = make_at2con(natoms, idef.il, idef.iparams, flexibleConstraintTreatment(bDynamics)); + const ListOfLists at2con = + make_at2con(natoms, idef.il, idef.iparams, flexibleConstraintTreatment(bDynamics)); const int ncon_tot = idef.il[F_CONSTR].nr / 3; @@ -2018,21 +2000,21 @@ void set_lincs(const t_idef& idef, const t_mdatoms& md, bool bDynamics, const t_ /* Skip the flexible constraints when not doing dynamics */ if (bDynamics || lenA != 0 || lenB != 0) { - assign_constraint(li, con, a1, a2, lenA, lenB, &at2con); + assign_constraint(li, con, a1, a2, lenA, lenB, at2con); if (li->ntask > 1 && !li->bTaskDep) { /* We can generate independent tasks. Check if we * need to assign connected constraints to our task. */ - check_assign_connected(li, iatom, idef, bDynamics, a1, a2, &at2con); + check_assign_connected(li, iatom, idef, bDynamics, a1, a2, at2con); } if (li->ntask > 1 && li->ncg_triangle > 0) { /* Ensure constraints in one triangle are assigned * to the same task. */ - check_assign_triangle(li, iatom, idef, bDynamics, con, a1, a2, &at2con); + check_assign_triangle(li, iatom, idef, bDynamics, con, a1, a2, at2con); } } } @@ -2096,13 +2078,11 @@ void set_lincs(const t_idef& idef, const t_mdatoms& md, bool bDynamics, const t_ li_task.tri_bits.resize(li_task.b1 - li_task.b0); } - set_matrix_indices(li, li_task, &at2con, bSortMatrix); + set_matrix_indices(li, li_task, at2con, bSortMatrix); } GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR } - done_blocka(&at2con); - if (cr->dd == nullptr) { /* Since the matrix is static, we should free some memory */ diff --git a/src/gromacs/mdlib/lincs.h b/src/gromacs/mdlib/lincs.h index dfc847127e..ef43c3a586 100644 --- a/src/gromacs/mdlib/lincs.h +++ b/src/gromacs/mdlib/lincs.h @@ -53,7 +53,6 @@ struct gmx_mtop_t; struct gmx_multisim_t; -struct t_blocka; struct t_commrec; struct t_idef; struct t_inputrec; @@ -65,9 +64,10 @@ namespace gmx { enum class ConstraintVariable : int; - -/* Abstract type for LINCS that is defined only in the file that uses it */ class Lincs; +template +class ListOfLists; + /*! \brief Return the data for determining constraint RMS relative deviations. */ ArrayRef lincs_rmsdData(Lincs* lincsd); @@ -76,13 +76,13 @@ ArrayRef lincs_rmsdData(Lincs* lincsd); real lincs_rmsd(const Lincs* lincsd); /*! \brief Initializes and returns the lincs data struct. */ -Lincs* init_lincs(FILE* fplog, - const gmx_mtop_t& mtop, - int nflexcon_global, - ArrayRef at2con, - bool bPLINCS, - int nIter, - int nProjOrder); +Lincs* init_lincs(FILE* fplog, + const gmx_mtop_t& mtop, + int nflexcon_global, + ArrayRef> atomsToConstraintsPerMolType, + bool bPLINCS, + int nIter, + int nProjOrder); /*! \brief Destructs the lincs object when it is not nullptr. */ void done_lincs(Lincs* li); diff --git a/src/gromacs/mdlib/tests/constrtestrunners.cpp b/src/gromacs/mdlib/tests/constrtestrunners.cpp index 710e9ac48c..ffd8f99e68 100644 --- a/src/gromacs/mdlib/tests/constrtestrunners.cpp +++ b/src/gromacs/mdlib/tests/constrtestrunners.cpp @@ -67,10 +67,10 @@ #include "gromacs/mdtypes/inputrec.h" #include "gromacs/mdtypes/mdatom.h" #include "gromacs/pbcutil/pbc.h" -#include "gromacs/topology/block.h" #include "gromacs/topology/idef.h" #include "gromacs/topology/ifunc.h" #include "gromacs/topology/topology.h" +#include "gromacs/utility/listoflists.h" #include "gromacs/utility/unique_cptr.h" #include "testutils/testasserts.h" @@ -115,7 +115,7 @@ void applyLincs(ConstraintsTestData* testData, t_pbc pbc) gmx_omp_nthreads_set(emntLINCS, 1); // Make blocka structure for faster LINCS setup - std::vector at2con_mt; + std::vector> at2con_mt; at2con_mt.reserve(testData->mtop_.moltype.size()); for (const gmx_moltype_t& moltype : testData->mtop_.moltype) { @@ -138,11 +138,6 @@ void applyLincs(ConstraintsTestData* testData, t_pbc pbc) &testData->nrnb_, maxwarn, &warncount_lincs); EXPECT_TRUE(success) << "Test failed with a false return value in LINCS."; EXPECT_EQ(warncount_lincs, 0) << "There were warnings in LINCS."; - for (auto& moltype : at2con_mt) - { - sfree(moltype.index); - sfree(moltype.a); - } done_lincs(lincsd); } diff --git a/src/gromacs/mdlib/updategroups.cpp b/src/gromacs/mdlib/updategroups.cpp index 5ee169a473..3aa3feb359 100644 --- a/src/gromacs/mdlib/updategroups.cpp +++ b/src/gromacs/mdlib/updategroups.cpp @@ -52,10 +52,10 @@ #include "gromacs/math/units.h" #include "gromacs/mdlib/constr.h" #include "gromacs/mdtypes/inputrec.h" -#include "gromacs/topology/block.h" #include "gromacs/topology/idef.h" #include "gromacs/topology/ifunc.h" #include "gromacs/topology/topology.h" +#include "gromacs/utility/listoflists.h" namespace gmx { @@ -196,15 +196,17 @@ static AtomIndexExtremes vsiteConstructRange(int a, const gmx_moltype_t& moltype } /*! \brief Returns the range of atoms constrained to atom \p a (including \p a itself) */ -static AtomIndexExtremes constraintAtomRange(int a, const t_blocka& at2con, const InteractionList& ilistConstraints) +static AtomIndexExtremes constraintAtomRange(int a, + const ListOfLists& at2con, + const InteractionList& ilistConstraints) { AtomIndexExtremes extremes = { a, a }; - for (int i = at2con.index[a]; i < at2con.index[a + 1]; i++) + for (const int constraint : at2con[a]) { for (int j = 0; j < 2; j++) { - int atomJ = ilistConstraints.iatoms[at2con.a[i] * 3 + 1 + j]; + int atomJ = ilistConstraints.iatoms[constraint * 3 + 1 + j]; extremes.minAtom = std::min(extremes.minAtom, atomJ); extremes.maxAtom = std::max(extremes.maxAtom, atomJ); } @@ -231,10 +233,10 @@ static std::vector buildIsParticleVsite(const gmx_moltype_t& moltype) } /*! \brief Returns the size of the update group starting at \p firstAtom or 0 when criteria (see updategroups.h) are not met */ -static int detectGroup(int firstAtom, - const gmx_moltype_t& moltype, - const t_blocka& at2con, - const InteractionList& ilistConstraints) +static int detectGroup(int firstAtom, + const gmx_moltype_t& moltype, + const ListOfLists& at2con, + const InteractionList& ilistConstraints) { /* We should be using moltype.atoms.atom[].ptype for checking whether * a particle is a vsite. But the test code can't fill t_atoms, @@ -243,7 +245,7 @@ static int detectGroup(int firstAtom, std::vector isParticleVsite = buildIsParticleVsite(moltype); /* A non-vsite atom without constraints is an update group by itself */ - if (!isParticleVsite[firstAtom] && at2con.index[firstAtom + 1] - at2con.index[firstAtom] == 0) + if (!isParticleVsite[firstAtom] && at2con[firstAtom].empty()) { return 1; } @@ -274,7 +276,7 @@ static int detectGroup(int firstAtom, } else { - int numConstraints = at2con.index[a + 1] - at2con.index[a]; + const int numConstraints = at2con[a].ssize(); if (numConstraints == 0) { /* We can not have unconstrained atoms in an update group */ @@ -356,8 +358,8 @@ static RangePartitioning makeUpdateGroups(const gmx_moltype_t& moltype, gmx::Arr ilistsCombined[F_CONSTR].iatoms = constraintsCombined.iatoms.data(); ilistsCombined[F_CONSTRNC].nr = 0; /* We "include" flexible constraints, but none are present (checked above) */ - t_blocka at2con = make_at2con(moltype.atoms.nr, ilistsCombined, iparams.data(), - FlexibleConstraintTreatment::Include); + const ListOfLists at2con = make_at2con(moltype.atoms.nr, ilistsCombined, iparams.data(), + FlexibleConstraintTreatment::Include); bool satisfiesCriteria = true; @@ -383,8 +385,6 @@ static RangePartitioning makeUpdateGroups(const gmx_moltype_t& moltype, gmx::Arr groups.clear(); } - done_blocka(&at2con); - return groups; } @@ -441,19 +441,19 @@ template static real constraintGroupRadius(const gmx_moltype_t& moltype, gmx::ArrayRef iparams, const int centralAtom, - const t_blocka& at2con, + const ListOfLists& at2con, const std::unordered_multimap& angleIndices, const real constraintLength, const real temperature) { - const int numConstraints = at2con.index[centralAtom + 1] - at2con.index[centralAtom]; + const int numConstraints = at2con[centralAtom].ssize(); GMX_RELEASE_ASSERT(numConstraints == numPartnerAtoms, "We expect as many constraints as partner atoms here"); std::array partnerAtoms; for (int i = 0; i < numPartnerAtoms; i++) { - const int ind = at2con.a[at2con.index[centralAtom] + i] * 3; + const int ind = at2con[centralAtom][i] * 3; if (ind >= moltype.ilist[F_CONSTR].size()) { /* This is a flexible constraint, we don't optimize for that */ @@ -597,7 +597,7 @@ static real computeMaxUpdateGroupRadius(const gmx_moltype_t& moltype, const InteractionList& settles = moltype.ilist[F_SETTLE]; - t_blocka at2con = make_at2con(moltype, iparams, FlexibleConstraintTreatment::Include); + const ListOfLists at2con = make_at2con(moltype, iparams, FlexibleConstraintTreatment::Include); const auto angleIndices = getAngleIndices(moltype); @@ -615,7 +615,7 @@ static real computeMaxUpdateGroupRadius(const gmx_moltype_t& moltype, int maxAtom = -1; for (int a : updateGroups.block(group)) { - int numConstraints = at2con.index[a + 1] - at2con.index[a]; + const int numConstraints = at2con[a].ssize(); if (numConstraints > maxNumConstraints) { maxNumConstraints = numConstraints; @@ -633,9 +633,10 @@ static real computeMaxUpdateGroupRadius(const gmx_moltype_t& moltype, int constraintType = -1; real maxConstraintLength = 0; real sumConstraintLengths = 0; - for (int i = at2con.index[maxAtom]; i < at2con.index[maxAtom + 1]; i++) + bool isFirstConstraint = true; + for (const int constraint : at2con[maxAtom]) { - int conIndex = at2con.a[i] * (1 + NRAL(F_CONSTR)); + int conIndex = constraint * (1 + NRAL(F_CONSTR)); int iparamsIndex; if (conIndex < moltype.ilist[F_CONSTR].size()) { @@ -646,9 +647,10 @@ static real computeMaxUpdateGroupRadius(const gmx_moltype_t& moltype, iparamsIndex = moltype.ilist[F_CONSTRNC].iatoms[conIndex - moltype.ilist[F_CONSTR].size()]; } - if (i == at2con.index[maxAtom]) + if (isFirstConstraint) { - constraintType = iparamsIndex; + constraintType = iparamsIndex; + isFirstConstraint = false; } else if (iparamsIndex != constraintType) { @@ -664,7 +666,7 @@ static real computeMaxUpdateGroupRadius(const gmx_moltype_t& moltype, sumConstraintLengths += constraintLength; } - int numConstraints = at2con.index[maxAtom + 1] - at2con.index[maxAtom]; + int numConstraints = at2con[maxAtom].ssize(); real radius; if (numConstraints == 1) { @@ -721,8 +723,6 @@ static real computeMaxUpdateGroupRadius(const gmx_moltype_t& moltype, maxRadius = std::max(maxRadius, dCAny); } - done_blocka(&at2con); - return maxRadius; } -- 2.22.0