From b749b01d1b9c7224bc7e2d50d1e2799586d98fcf Mon Sep 17 00:00:00 2001 From: Mark Abraham Date: Mon, 31 May 2021 11:05:50 +0200 Subject: [PATCH] Split localtopology from domdec_topology This is a separable responsibility and can go in its own translation unit. Refs #3887 --- src/gromacs/domdec/domdec.h | 17 - src/gromacs/domdec/domdec_topology.cpp | 959 ---------------------- src/gromacs/domdec/localtopology.cpp | 1031 ++++++++++++++++++++++++ src/gromacs/domdec/localtopology.h | 79 ++ src/gromacs/domdec/mdsetup.cpp | 1 + src/gromacs/domdec/partition.cpp | 2 +- 6 files changed, 1112 insertions(+), 977 deletions(-) create mode 100644 src/gromacs/domdec/localtopology.cpp create mode 100644 src/gromacs/domdec/localtopology.h diff --git a/src/gromacs/domdec/domdec.h b/src/gromacs/domdec/domdec.h index 55644322b7..e6564bf752 100644 --- a/src/gromacs/domdec/domdec.h +++ b/src/gromacs/domdec/domdec.h @@ -241,23 +241,6 @@ gmx::ArrayRef dd_constraints_nlocalatoms(const gmx_domdec_t* dd); /* In domdec_top.c */ -/*! \brief Generate the local topology and virtual site data - * - * \returns Total count of bonded interactions in the local topology on this domain */ -int dd_make_local_top(struct gmx_domdec_t* dd, - struct gmx_domdec_zones_t* zones, - int npbcdim, - matrix box, - rvec cellsize_min, - const ivec npulse, - t_forcerec* fr, - gmx::ArrayRef coordinates, - const gmx_mtop_t& top, - gmx_localtop_t* ltop); - -/*! \brief Sort ltop->ilist when we are doing free energy. */ -void dd_sort_local_top(const gmx_domdec_t& dd, const t_mdatoms* mdatoms, gmx_localtop_t* ltop); - /*! \brief Construct local state */ void dd_init_local_state(const gmx_domdec_t& dd, const t_state* state_global, t_state* local_state); diff --git a/src/gromacs/domdec/domdec_topology.cpp b/src/gromacs/domdec/domdec_topology.cpp index 2586694f77..289a059871 100644 --- a/src/gromacs/domdec/domdec_topology.cpp +++ b/src/gromacs/domdec/domdec_topology.cpp @@ -94,970 +94,11 @@ using gmx::ArrayRef; using gmx::DDBondedChecking; -using gmx::ListOfLists; using gmx::RVec; /*! \brief The number of integer item in the local state, used for broadcasting of the state */ #define NITEM_DD_INIT_LOCAL_STATE 5 -/*! \brief Store a vsite interaction at the end of \p il - * - * This routine is very similar to add_ifunc, but vsites interactions - * have more work to do than other kinds of interactions, and the - * complex way nral (and thus vector contents) depends on ftype - * confuses static analysis tools unless we fuse the vsite - * atom-indexing organization code with the ifunc-adding code, so that - * they can see that nral is the same value. */ -static inline void add_ifunc_for_vsites(t_iatom* tiatoms, - const gmx_ga2la_t& ga2la, - int nral, - gmx_bool bHomeA, - int a, - int a_gl, - int a_mol, - const t_iatom* iatoms, - InteractionList* il) -{ - /* Copy the type */ - tiatoms[0] = iatoms[0]; - - if (bHomeA) - { - /* We know the local index of the first atom */ - tiatoms[1] = a; - } - else - { - /* Convert later in make_local_vsites */ - tiatoms[1] = -a_gl - 1; - } - - GMX_ASSERT(nral >= 2 && nral <= 5, "Invalid nral for vsites"); - for (int k = 2; k < 1 + nral; k++) - { - int ak_gl = a_gl + iatoms[k] - a_mol; - if (const int* homeIndex = ga2la.findHome(ak_gl)) - { - tiatoms[k] = *homeIndex; - } - else - { - /* Copy the global index, convert later in make_local_vsites */ - tiatoms[k] = -(ak_gl + 1); - } - // Note that ga2la_get_home always sets the third parameter if - // it returns TRUE - } - 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, - gmx::ArrayRef iatoms, - const t_iparams* ip_in, - InteractionDefinitions* idef) -{ - /* This position restraint has not been added yet, - * so it's index is the current number of position restraints. - */ - const int n = idef->il[F_POSRES].size() / 2; - - /* Get the position restraint coordinates from the molblock */ - 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]]; - 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]; - } - 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]; - } - /* 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, - gmx::ArrayRef iatoms, - const t_iparams* ip_in, - InteractionDefinitions* idef) -{ - /* This flat-bottom position restraint has not been added yet, - * so it's index is the current number of position restraints. - */ - const int n = idef->il[F_FBPOSRES].size() / 2; - - /* Get the position restraint coordinats from the molblock */ - 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]; - - /* Note: no B-type for flat-bottom posres */ - - /* 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 */ -static void add_vsite(const gmx_ga2la_t& ga2la, - gmx::ArrayRef index, - gmx::ArrayRef rtil, - int ftype, - int nral, - gmx_bool bHomeA, - int a, - int a_gl, - int a_mol, - const t_iatom* iatoms, - InteractionDefinitions* idef) -{ - t_iatom tiatoms[1 + MAXATOMLIST]; - - /* Add this interaction to the local topology */ - add_ifunc_for_vsites(tiatoms, ga2la, nral, bHomeA, a, a_gl, a_mol, iatoms, &idef->il[ftype]); - - if (iatoms[1 + nral]) - { - /* Check for recursion */ - for (int k = 2; k < 1 + nral; k++) - { - if ((iatoms[1 + nral] & (2 << k)) && (tiatoms[k] < 0)) - { - /* This construction atoms is a vsite and not a home atom */ - if (gmx_debug_at) - { - fprintf(debug, - "Constructing atom %d of vsite atom %d is a vsite and non-home\n", - iatoms[k] + 1, - a_mol + 1); - } - /* Find the vsite construction */ - - /* Check all interactions assigned to this atom */ - int j = index[iatoms[k]]; - while (j < index[iatoms[k] + 1]) - { - int ftype_r = rtil[j++]; - int nral_r = NRAL(ftype_r); - if (interaction_function[ftype_r].flags & IF_VSITE) - { - /* Add this vsite (recursion) */ - add_vsite(ga2la, - index, - rtil, - ftype_r, - nral_r, - FALSE, - -1, - a_gl + iatoms[k] - iatoms[1], - iatoms[k], - rtil.data() + j, - idef); - } - j += 1 + nral_rt(ftype_r); - } - } - } - } -} - -/*! \brief Returns the squared distance between atoms \p i and \p j */ -static real dd_dist2(const t_pbc* pbc_null, ArrayRef coordinates, const int i, const int j) -{ - rvec dx; - - if (pbc_null) - { - pbc_dx_aiuc(pbc_null, coordinates[i], coordinates[j], dx); - } - else - { - rvec_sub(coordinates[i], coordinates[j], dx); - } - - return norm2(dx); -} - -/*! \brief Append t_idef structures 1 to nsrc in src to *dest */ -static void combine_idef(InteractionDefinitions* dest, gmx::ArrayRef src) -{ - for (int ftype = 0; ftype < F_NRE; ftype++) - { - int n = 0; - for (gmx::index s = 1; s < src.ssize(); s++) - { - n += src[s].idef.il[ftype].size(); - } - if (n > 0) - { - for (gmx::index s = 1; s < src.ssize(); s++) - { - 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].size() / 2; - std::vector& iparams_dest = - (ftype == F_POSRES ? dest->iparams_posres : dest->iparams_fbposres); - - /* 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].size() / 2; - } - - for (gmx::index s = 1; s < src.ssize(); s++) - { - 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()); - - /* 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; - nposres++; - } - } - GMX_RELEASE_ASSERT( - int(iparams_dest.size()) == nposres, - "The number of parameters should match the number of restraints"); - } - } - } -} - -//! Options for assigning interactions for atoms -enum class InteractionConnectivity -{ - Intramolecular, //!< Only intra-molecular interactions - Intermolecular //!< Only inter-molecular interactions -}; - -/*! \brief Determine whether the local domain has responsibility for - * any of the bonded interactions for local atom \p atomIndex - * and assign those to the local domain. - * - * \returns The total number of bonded interactions for this atom for - * which this domain is responsible. - */ -template -static inline int assignInteractionsForAtom(const int atomIndex, - const int globalAtomIndex, - const int atomIndexInMolecule, - gmx::ArrayRef index, - gmx::ArrayRef rtil, - const int ind_start, - const int ind_end, - const gmx_ga2la_t& ga2la, - const gmx_domdec_zones_t* zones, - const bool checkDistanceMultiBody, - const ivec rcheck, - const bool checkDistanceTwoBody, - const real cutoffSquared, - const t_pbc* pbc_null, - ArrayRef coordinates, - InteractionDefinitions* idef, - const int iz, - const DDBondedChecking ddBondedChecking) -{ - gmx::ArrayRef iZones = zones->iZones; - - int numBondedInteractions = 0; - - int j = ind_start; - while (j < ind_end) - { - t_iatom tiatoms[1 + MAXATOMLIST]; - - const int ftype = rtil[j++]; - auto iatoms = gmx::constArrayRefFromArray(rtil.data() + j, rtil.size() - j); - const int nral = NRAL(ftype); - if (interaction_function[ftype].flags & IF_VSITE) - { - GMX_ASSERT(interactionConnectivity == InteractionConnectivity::Intramolecular, - "All vsites should be intramolecular"); - - /* The vsite construction goes where the vsite itself is */ - if (iz == 0) - { - add_vsite(ga2la, - index, - rtil, - ftype, - nral, - TRUE, - atomIndex, - globalAtomIndex, - atomIndexInMolecule, - iatoms.data(), - idef); - } - } - else - { - bool bUse = false; - - /* Copy the type */ - tiatoms[0] = iatoms[0]; - - if (nral == 1) - { - GMX_ASSERT(interactionConnectivity == InteractionConnectivity::Intramolecular, - "All interactions that involve a single atom are intramolecular"); - - /* Assign single-body interactions to the home zone. - * Position restraints are not handled here, but separately. - */ - if (iz == 0 && !(ftype == F_POSRES || ftype == F_FBPOSRES)) - { - bUse = true; - tiatoms[1] = atomIndex; - } - } - else if (nral == 2) - { - /* This is a two-body interaction, we can assign - * analogous to the non-bonded assignments. - */ - const int k_gl = (interactionConnectivity == InteractionConnectivity::Intramolecular) - ? - /* Get the global index using the offset in the molecule */ - (globalAtomIndex + iatoms[2] - atomIndexInMolecule) - : iatoms[2]; - if (const auto* entry = ga2la.find(k_gl)) - { - int kz = entry->cell; - if (kz >= zones->n) - { - kz -= zones->n; - } - /* Check zone interaction assignments */ - bUse = ((iz < iZones.ssize() && iz <= kz && iZones[iz].jZoneRange.isInRange(kz)) - || (kz < iZones.ssize() && iz > kz && iZones[kz].jZoneRange.isInRange(iz))); - if (bUse) - { - GMX_ASSERT(ftype != F_CONSTR || (iz == 0 && kz == 0), - "Constraint assigned here should only involve home atoms"); - - tiatoms[1] = atomIndex; - tiatoms[2] = entry->la; - /* If necessary check the cgcm distance */ - if (checkDistanceTwoBody - && dd_dist2(pbc_null, coordinates, tiatoms[1], tiatoms[2]) >= cutoffSquared) - { - bUse = false; - } - } - } - else - { - bUse = false; - } - } - else - { - /* Assign this multi-body bonded interaction to - * the local domain if we have all the atoms involved - * (local or communicated) and the minimum zone shift - * in each dimension is zero, for dimensions - * with 2 DD cells an extra check may be necessary. - */ - ivec k_zero, k_plus; - bUse = true; - clear_ivec(k_zero); - clear_ivec(k_plus); - for (int k = 1; k <= nral && bUse; k++) - { - const int k_gl = - (interactionConnectivity == InteractionConnectivity::Intramolecular) - ? - /* Get the global index using the offset in the molecule */ - (globalAtomIndex + iatoms[k] - atomIndexInMolecule) - : iatoms[k]; - const auto* entry = ga2la.find(k_gl); - if (entry == nullptr || entry->cell >= zones->n) - { - /* We do not have this atom of this interaction - * locally, or it comes from more than one cell - * away. - */ - bUse = FALSE; - } - else - { - tiatoms[k] = entry->la; - for (int d = 0; d < DIM; d++) - { - if (zones->shift[entry->cell][d] == 0) - { - k_zero[d] = k; - } - else - { - k_plus[d] = k; - } - } - } - } - bUse = (bUse && (k_zero[XX] != 0) && (k_zero[YY] != 0) && (k_zero[ZZ] != 0)); - if (checkDistanceMultiBody) - { - for (int d = 0; (d < DIM && bUse); d++) - { - /* Check if the distance falls within - * the cut-off to avoid possible multiple - * assignments of bonded interactions. - */ - if (rcheck[d] && k_plus[d] - && dd_dist2(pbc_null, coordinates, tiatoms[k_zero[d]], tiatoms[k_plus[d]]) - >= cutoffSquared) - { - bUse = FALSE; - } - } - } - } - if (bUse) - { - /* Add this interaction to the local topology */ - idef->il[ftype].push_back(tiatoms[0], nral, tiatoms + 1); - /* Sum so we can check in global_stat - * if we have everything. - */ - if (ddBondedChecking == DDBondedChecking::All - || !(interaction_function[ftype].flags & IF_LIMZERO)) - { - numBondedInteractions++; - } - } - } - j += 1 + nral_rt(ftype); - } - - return numBondedInteractions; -} - -/*! \brief Determine whether the local domain has responsibility for - * any of the position restraints for local atom \p atomIndex - * and assign those to the local domain. - * - * \returns The total number of bonded interactions for this atom for - * which this domain is responsible. - */ -static inline int assignPositionRestraintsForAtom(const int atomIndex, - const int mol, - const int atomIndexInMolecule, - const int numAtomsInMolecule, - gmx::ArrayRef rtil, - const int ind_start, - const int ind_end, - const gmx_molblock_t& molb, - const t_iparams* ip_in, - InteractionDefinitions* idef) -{ - constexpr int nral = 1; - - int numBondedInteractions = 0; - - int j = ind_start; - while (j < ind_end) - { - const int ftype = rtil[j++]; - auto iatoms = gmx::constArrayRefFromArray(rtil.data() + j, rtil.size() - j); - - if (ftype == F_POSRES || ftype == F_FBPOSRES) - { - std::array tiatoms = { iatoms[0], atomIndex }; - if (ftype == F_POSRES) - { - add_posres(mol, atomIndexInMolecule, numAtomsInMolecule, molb, tiatoms, ip_in, idef); - } - else - { - add_fbposres(mol, atomIndexInMolecule, numAtomsInMolecule, molb, tiatoms, ip_in, idef); - } - idef->il[ftype].push_back(tiatoms[0], nral, tiatoms.data() + 1); - numBondedInteractions++; - } - j += 1 + nral_rt(ftype); - } - - return numBondedInteractions; -} - -/*! \brief This function looks up and assigns bonded interactions for zone iz. - * - * With thread parallelizing each thread acts on a different atom range: - * at_start to at_end. - */ -static int make_bondeds_zone(gmx_reverse_top_t* rt, - ArrayRef globalAtomIndices, - const gmx_ga2la_t& ga2la, - const gmx_domdec_zones_t* zones, - const std::vector& molb, - const bool checkDistanceMultiBody, - const ivec rcheck, - const bool checkDistanceTwoBody, - const real cutoffSquared, - const t_pbc* pbc_null, - ArrayRef coordinates, - const t_iparams* ip_in, - InteractionDefinitions* idef, - int izone, - const gmx::Range& atomRange) -{ - int mb = 0; - int mt = 0; - int mol = 0; - int atomIndexInMolecule = 0; - - const auto ddBondedChecking = rt->options().ddBondedChecking; - - int numBondedInteractions = 0; - - for (int i : atomRange) - { - /* Get the global atom number */ - const int globalAtomIndex = globalAtomIndices[i]; - global_atomnr_to_moltype_ind( - rt->molblockIndices(), globalAtomIndex, &mb, &mt, &mol, &atomIndexInMolecule); - /* Check all intramolecular interactions assigned to this atom */ - gmx::ArrayRef index = rt->interactionListForMoleculeType(mt).index; - gmx::ArrayRef rtil = rt->interactionListForMoleculeType(mt).il; - - numBondedInteractions += assignInteractionsForAtom( - i, - globalAtomIndex, - atomIndexInMolecule, - index, - rtil, - index[atomIndexInMolecule], - index[atomIndexInMolecule + 1], - ga2la, - zones, - checkDistanceMultiBody, - rcheck, - checkDistanceTwoBody, - cutoffSquared, - pbc_null, - coordinates, - idef, - izone, - ddBondedChecking); - - // Assign position restraints, when present, for the home zone - if (izone == 0 && rt->hasPositionRestraints()) - { - numBondedInteractions += assignPositionRestraintsForAtom( - i, - mol, - atomIndexInMolecule, - rt->interactionListForMoleculeType(mt).numAtomsInMolecule, - rtil, - index[atomIndexInMolecule], - index[atomIndexInMolecule + 1], - molb[mb], - ip_in, - idef); - } - - if (rt->hasIntermolecularInteractions()) - { - /* Check all intermolecular interactions assigned to this atom */ - index = rt->interactionListForIntermolecularInteractions().index; - rtil = rt->interactionListForIntermolecularInteractions().il; - - numBondedInteractions += assignInteractionsForAtom( - i, - -1, // not used - -1, // not used - index, - rtil, - index[globalAtomIndex], - index[globalAtomIndex + 1], - ga2la, - zones, - checkDistanceMultiBody, - rcheck, - checkDistanceTwoBody, - cutoffSquared, - pbc_null, - coordinates, - idef, - izone, - ddBondedChecking); - } - } - - return numBondedInteractions; -} - -/*! \brief Set the exclusion data for i-zone \p iz */ -static void make_exclusions_zone(ArrayRef globalAtomIndices, - const gmx_ga2la_t& ga2la, - gmx_domdec_zones_t* zones, - ArrayRef molblockIndices, - const std::vector& moltype, - const int* atomInfo, - ListOfLists* lexcls, - int iz, - int at_start, - int at_end, - const gmx::ArrayRef intermolecularExclusionGroup) -{ - const auto& jAtomRange = zones->iZones[iz].jAtomRange; - - const gmx::index oldNumLists = lexcls->ssize(); - - std::vector exclusionsForAtom; - for (int at = at_start; at < at_end; at++) - { - exclusionsForAtom.clear(); - - if (GET_CGINFO_EXCL_INTER(atomInfo[at])) - { - int mb = 0; - int mt = 0; - int mol = 0; - int a_mol = 0; - - /* Copy the exclusions from the global top */ - int a_gl = globalAtomIndices[at]; - global_atomnr_to_moltype_ind(molblockIndices, a_gl, &mb, &mt, &mol, &a_mol); - const auto excls = moltype[mt].excls[a_mol]; - for (const int aj_mol : excls) - { - if (const auto* jEntry = ga2la.find(a_gl + aj_mol - a_mol)) - { - /* This check is not necessary, but it can reduce - * the number of exclusions in the list, which in turn - * can speed up the pair list construction a bit. - */ - if (jAtomRange.isInRange(jEntry->la)) - { - exclusionsForAtom.push_back(jEntry->la); - } - } - } - } - - bool isExcludedAtom = !intermolecularExclusionGroup.empty() - && std::find(intermolecularExclusionGroup.begin(), - intermolecularExclusionGroup.end(), - globalAtomIndices[at]) - != intermolecularExclusionGroup.end(); - - if (isExcludedAtom) - { - for (int qmAtomGlobalIndex : intermolecularExclusionGroup) - { - if (const auto* entry = ga2la.find(qmAtomGlobalIndex)) - { - exclusionsForAtom.push_back(entry->la); - } - } - } - - /* Append the exclusions for this atom to the topology */ - lexcls->pushBack(exclusionsForAtom); - } - - GMX_RELEASE_ASSERT( - lexcls->ssize() - oldNumLists == at_end - at_start, - "The number of exclusion list should match the number of atoms in the range"); -} - -/*! \brief Generate and store all required local bonded interactions in \p idef and local exclusions in \p lexcls - * - * \returns Total count of bonded interactions in the local topology on this domain */ -static int make_local_bondeds_excls(gmx_domdec_t* dd, - gmx_domdec_zones_t* zones, - const gmx_mtop_t& mtop, - const int* atomInfo, - const bool checkDistanceMultiBody, - const ivec rcheck, - const gmx_bool checkDistanceTwoBody, - const real cutoff, - const t_pbc* pbc_null, - ArrayRef coordinates, - InteractionDefinitions* idef, - ListOfLists* lexcls) -{ - int nzone_bondeds = 0; - - if (dd->reverse_top->hasInterAtomicInteractions()) - { - nzone_bondeds = zones->n; - } - else - { - /* Only single charge group (or atom) molecules, so interactions don't - * cross zone boundaries and we only need to assign in the home zone. - */ - nzone_bondeds = 1; - } - - /* We only use exclusions from i-zones to i- and j-zones */ - const int numIZonesForExclusions = (dd->haveExclusions ? zones->iZones.size() : 0); - - gmx_reverse_top_t* rt = dd->reverse_top.get(); - - const real cutoffSquared = gmx::square(cutoff); - - /* Clear the counts */ - idef->clear(); - int numBondedInteractions = 0; - - lexcls->clear(); - - for (int izone = 0; izone < nzone_bondeds; izone++) - { - const int cg0 = zones->cg_range[izone]; - const int cg1 = zones->cg_range[izone + 1]; - - gmx::ArrayRef threadWorkObjects = rt->threadWorkObjects(); - const int numThreads = threadWorkObjects.size(); -#pragma omp parallel for num_threads(numThreads) schedule(static) - for (int thread = 0; thread < numThreads; thread++) - { - try - { - InteractionDefinitions* idef_t = nullptr; - - int cg0t = cg0 + ((cg1 - cg0) * thread) / numThreads; - int cg1t = cg0 + ((cg1 - cg0) * (thread + 1)) / numThreads; - - if (thread == 0) - { - idef_t = idef; - } - else - { - idef_t = &threadWorkObjects[thread].idef; - idef_t->clear(); - } - - threadWorkObjects[thread].numBondedInteractions = - make_bondeds_zone(rt, - dd->globalAtomIndices, - *dd->ga2la, - zones, - mtop.molblock, - checkDistanceMultiBody, - rcheck, - checkDistanceTwoBody, - cutoffSquared, - pbc_null, - coordinates, - idef->iparams.data(), - idef_t, - izone, - gmx::Range(cg0t, cg1t)); - - if (izone < numIZonesForExclusions) - { - ListOfLists* excl_t = nullptr; - if (thread == 0) - { - // Thread 0 stores exclusions directly in the final storage - excl_t = lexcls; - } - else - { - // Threads > 0 store in temporary storage, starting at list index 0 - excl_t = &threadWorkObjects[thread].excl; - excl_t->clear(); - } - - /* No charge groups and no distance check required */ - make_exclusions_zone(dd->globalAtomIndices, - *dd->ga2la, - zones, - rt->molblockIndices(), - mtop.moltype, - atomInfo, - excl_t, - izone, - cg0t, - cg1t, - mtop.intermolecularExclusionGroup); - } - } - GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR - } - - if (threadWorkObjects.size() > 1) - { - combine_idef(idef, threadWorkObjects); - } - - for (const thread_work_t& th_work : threadWorkObjects) - { - numBondedInteractions += th_work.numBondedInteractions; - } - - if (izone < numIZonesForExclusions) - { - for (std::size_t th = 1; th < threadWorkObjects.size(); th++) - { - lexcls->appendListOfLists(threadWorkObjects[th].excl); - } - } - } - - if (debug) - { - fprintf(debug, "We have %d exclusions\n", lexcls->numElements()); - } - - return numBondedInteractions; -} - -int dd_make_local_top(gmx_domdec_t* dd, - gmx_domdec_zones_t* zones, - int npbcdim, - matrix box, - rvec cellsize_min, - const ivec npulse, - t_forcerec* fr, - ArrayRef coordinates, - const gmx_mtop_t& mtop, - gmx_localtop_t* ltop) -{ - real rc = -1; - ivec rcheck; - t_pbc pbc, *pbc_null = nullptr; - - if (debug) - { - fprintf(debug, "Making local topology\n"); - } - - bool checkDistanceMultiBody = false; - bool checkDistanceTwoBody = false; - - if (dd->reverse_top->hasInterAtomicInteractions()) - { - /* We need to check to which cell bondeds should be assigned */ - rc = dd_cutoff_twobody(dd); - if (debug) - { - fprintf(debug, "Two-body bonded cut-off distance is %g\n", rc); - } - - /* Should we check distances when assigning bonded interactions? */ - for (int d = 0; d < DIM; d++) - { - rcheck[d] = FALSE; - /* Only need to check for dimensions where the part of the box - * that is not communicated is smaller than the cut-off. - */ - if (d < npbcdim && dd->numCells[d] > 1 - && (dd->numCells[d] - npulse[d]) * cellsize_min[d] < 2 * rc) - { - if (dd->numCells[d] == 2) - { - rcheck[d] = TRUE; - checkDistanceMultiBody = TRUE; - } - /* Check for interactions between two atoms, - * where we can allow interactions up to the cut-off, - * instead of up to the smallest cell dimension. - */ - checkDistanceTwoBody = TRUE; - } - if (debug) - { - fprintf(debug, - "dim %d cellmin %f bonded rcheck[%d] = %d, checkDistanceTwoBody = %s\n", - d, - cellsize_min[d], - d, - rcheck[d], - gmx::boolToString(checkDistanceTwoBody)); - } - } - if (checkDistanceMultiBody || checkDistanceTwoBody) - { - if (fr->bMolPBC) - { - pbc_null = set_pbc_dd(&pbc, fr->pbcType, dd->numCells, TRUE, box); - } - else - { - pbc_null = nullptr; - } - } - } - - int numBondedInteractionsToReduce = make_local_bondeds_excls(dd, - zones, - mtop, - fr->atomInfo.data(), - checkDistanceMultiBody, - rcheck, - checkDistanceTwoBody, - rc, - pbc_null, - coordinates, - <op->idef, - <op->excls); - - /* The ilist is not sorted yet, - * we can only do this when we have the charge arrays. - */ - ltop->idef.ilsort = ilsortUNKNOWN; - - return numBondedInteractionsToReduce; -} - -void dd_sort_local_top(const gmx_domdec_t& dd, const t_mdatoms* mdatoms, gmx_localtop_t* ltop) -{ - if (dd.reverse_top->doSorting()) - { - gmx_sort_ilist_fe(<op->idef, mdatoms->chargeA, mdatoms->chargeB); - } - else - { - ltop->idef.ilsort = ilsortNO_FE; - } -} - void dd_init_local_state(const gmx_domdec_t& dd, const t_state* state_global, t_state* state_local) { int buf[NITEM_DD_INIT_LOCAL_STATE]; diff --git a/src/gromacs/domdec/localtopology.cpp b/src/gromacs/domdec/localtopology.cpp new file mode 100644 index 0000000000..ace7b43a1b --- /dev/null +++ b/src/gromacs/domdec/localtopology.cpp @@ -0,0 +1,1031 @@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright (c) 2006 - 2014, The GROMACS development team. + * Copyright (c) 2015,2016,2017,2018,2019,2020,2021, 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. + * + * GROMACS is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + * + * GROMACS is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with GROMACS; if not, see + * http://www.gnu.org/licenses, or write to the Free Software Foundation, + * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + * + * If you want to redistribute modifications to GROMACS, please + * consider that scientific software is very special. Version + * control is crucial - bugs must be traceable. We will be happy to + * consider code for inclusion in the official distribution, but + * derived work must not be called official GROMACS. Details are found + * in the README & COPYING files - if they are missing, get the + * official version at http://www.gromacs.org. + * + * To help us fund GROMACS development, we humbly ask that you cite + * the research papers on the package. Check out http://www.gromacs.org. + */ + +/*! \internal \file + * + * \brief This file defines functions used by the domdec module while + * building topologies local to a domain + * + * \author Berk Hess + * \ingroup module_domdec + */ + +#include "gmxpre.h" + +#include "gromacs/domdec/localtopology.h" + +#include + +#include "gromacs/domdec/domdec_internal.h" +#include "gromacs/domdec/domdec_struct.h" +#include "gromacs/domdec/ga2la.h" +#include "gromacs/domdec/options.h" +#include "gromacs/domdec/reversetopology.h" +#include "gromacs/math/vec.h" +#include "gromacs/mdtypes/forcerec.h" +#include "gromacs/mdtypes/mdatom.h" +#include "gromacs/pbcutil/pbc.h" +#include "gromacs/topology/idef.h" +#include "gromacs/topology/mtop_util.h" +#include "gromacs/topology/topsort.h" +#include "gromacs/utility/arrayref.h" +#include "gromacs/utility/fatalerror.h" +#include "gromacs/utility/gmxassert.h" +#include "gromacs/utility/listoflists.h" +#include "gromacs/utility/strconvert.h" + +using gmx::ArrayRef; +using gmx::DDBondedChecking; +using gmx::ListOfLists; +using gmx::RVec; + +/*! \brief Store a vsite interaction at the end of \p il + * + * This routine is very similar to add_ifunc, but vsites interactions + * have more work to do than other kinds of interactions, and the + * complex way nral (and thus vector contents) depends on ftype + * confuses static analysis tools unless we fuse the vsite + * atom-indexing organization code with the ifunc-adding code, so that + * they can see that nral is the same value. */ +static inline void add_ifunc_for_vsites(t_iatom* tiatoms, + const gmx_ga2la_t& ga2la, + int nral, + gmx_bool bHomeA, + int a, + int a_gl, + int a_mol, + const t_iatom* iatoms, + InteractionList* il) +{ + /* Copy the type */ + tiatoms[0] = iatoms[0]; + + if (bHomeA) + { + /* We know the local index of the first atom */ + tiatoms[1] = a; + } + else + { + /* Convert later in make_local_vsites */ + tiatoms[1] = -a_gl - 1; + } + + GMX_ASSERT(nral >= 2 && nral <= 5, "Invalid nral for vsites"); + for (int k = 2; k < 1 + nral; k++) + { + int ak_gl = a_gl + iatoms[k] - a_mol; + if (const int* homeIndex = ga2la.findHome(ak_gl)) + { + tiatoms[k] = *homeIndex; + } + else + { + /* Copy the global index, convert later in make_local_vsites */ + tiatoms[k] = -(ak_gl + 1); + } + // Note that ga2la_get_home always sets the third parameter if + // it returns TRUE + } + 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, + gmx::ArrayRef iatoms, + const t_iparams* ip_in, + InteractionDefinitions* idef) +{ + /* This position restraint has not been added yet, + * so it's index is the current number of position restraints. + */ + const int n = idef->il[F_POSRES].size() / 2; + + /* Get the position restraint coordinates from the molblock */ + 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]]; + 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]; + } + 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]; + } + /* 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, + gmx::ArrayRef iatoms, + const t_iparams* ip_in, + InteractionDefinitions* idef) +{ + /* This flat-bottom position restraint has not been added yet, + * so it's index is the current number of position restraints. + */ + const int n = idef->il[F_FBPOSRES].size() / 2; + + /* Get the position restraint coordinats from the molblock */ + 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]; + + /* Note: no B-type for flat-bottom posres */ + + /* 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 */ +static void add_vsite(const gmx_ga2la_t& ga2la, + gmx::ArrayRef index, + gmx::ArrayRef rtil, + int ftype, + int nral, + gmx_bool bHomeA, + int a, + int a_gl, + int a_mol, + const t_iatom* iatoms, + InteractionDefinitions* idef) +{ + t_iatom tiatoms[1 + MAXATOMLIST]; + + /* Add this interaction to the local topology */ + add_ifunc_for_vsites(tiatoms, ga2la, nral, bHomeA, a, a_gl, a_mol, iatoms, &idef->il[ftype]); + + if (iatoms[1 + nral]) + { + /* Check for recursion */ + for (int k = 2; k < 1 + nral; k++) + { + if ((iatoms[1 + nral] & (2 << k)) && (tiatoms[k] < 0)) + { + /* This construction atoms is a vsite and not a home atom */ + if (gmx_debug_at) + { + fprintf(debug, + "Constructing atom %d of vsite atom %d is a vsite and non-home\n", + iatoms[k] + 1, + a_mol + 1); + } + /* Find the vsite construction */ + + /* Check all interactions assigned to this atom */ + int j = index[iatoms[k]]; + while (j < index[iatoms[k] + 1]) + { + int ftype_r = rtil[j++]; + int nral_r = NRAL(ftype_r); + if (interaction_function[ftype_r].flags & IF_VSITE) + { + /* Add this vsite (recursion) */ + add_vsite(ga2la, + index, + rtil, + ftype_r, + nral_r, + FALSE, + -1, + a_gl + iatoms[k] - iatoms[1], + iatoms[k], + rtil.data() + j, + idef); + } + j += 1 + nral_rt(ftype_r); + } + } + } + } +} + +/*! \brief Returns the squared distance between atoms \p i and \p j */ +static real dd_dist2(const t_pbc* pbc_null, ArrayRef coordinates, const int i, const int j) +{ + rvec dx; + + if (pbc_null) + { + pbc_dx_aiuc(pbc_null, coordinates[i], coordinates[j], dx); + } + else + { + rvec_sub(coordinates[i], coordinates[j], dx); + } + + return norm2(dx); +} + +/*! \brief Append t_idef structures 1 to nsrc in src to *dest */ +static void combine_idef(InteractionDefinitions* dest, gmx::ArrayRef src) +{ + for (int ftype = 0; ftype < F_NRE; ftype++) + { + int n = 0; + for (gmx::index s = 1; s < src.ssize(); s++) + { + n += src[s].idef.il[ftype].size(); + } + if (n > 0) + { + for (gmx::index s = 1; s < src.ssize(); s++) + { + 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].size() / 2; + std::vector& iparams_dest = + (ftype == F_POSRES ? dest->iparams_posres : dest->iparams_fbposres); + + /* 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].size() / 2; + } + + for (gmx::index s = 1; s < src.ssize(); s++) + { + 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()); + + /* 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; + nposres++; + } + } + GMX_RELEASE_ASSERT( + int(iparams_dest.size()) == nposres, + "The number of parameters should match the number of restraints"); + } + } + } +} + +//! Options for assigning interactions for atoms +enum class InteractionConnectivity +{ + Intramolecular, //!< Only intra-molecular interactions + Intermolecular //!< Only inter-molecular interactions +}; + +/*! \brief Determine whether the local domain has responsibility for + * any of the bonded interactions for local atom \p atomIndex + * and assign those to the local domain. + * + * \returns The total number of bonded interactions for this atom for + * which this domain is responsible. + */ +template +static inline int assignInteractionsForAtom(const int atomIndex, + const int globalAtomIndex, + const int atomIndexInMolecule, + gmx::ArrayRef index, + gmx::ArrayRef rtil, + const int ind_start, + const int ind_end, + const gmx_ga2la_t& ga2la, + const gmx_domdec_zones_t* zones, + const bool checkDistanceMultiBody, + const ivec rcheck, + const bool checkDistanceTwoBody, + const real cutoffSquared, + const t_pbc* pbc_null, + ArrayRef coordinates, + InteractionDefinitions* idef, + const int iz, + const DDBondedChecking ddBondedChecking) +{ + gmx::ArrayRef iZones = zones->iZones; + + int numBondedInteractions = 0; + + int j = ind_start; + while (j < ind_end) + { + t_iatom tiatoms[1 + MAXATOMLIST]; + + const int ftype = rtil[j++]; + auto iatoms = gmx::constArrayRefFromArray(rtil.data() + j, rtil.size() - j); + const int nral = NRAL(ftype); + if (interaction_function[ftype].flags & IF_VSITE) + { + GMX_ASSERT(interactionConnectivity == InteractionConnectivity::Intramolecular, + "All vsites should be intramolecular"); + + /* The vsite construction goes where the vsite itself is */ + if (iz == 0) + { + add_vsite(ga2la, + index, + rtil, + ftype, + nral, + TRUE, + atomIndex, + globalAtomIndex, + atomIndexInMolecule, + iatoms.data(), + idef); + } + } + else + { + bool bUse = false; + + /* Copy the type */ + tiatoms[0] = iatoms[0]; + + if (nral == 1) + { + GMX_ASSERT(interactionConnectivity == InteractionConnectivity::Intramolecular, + "All interactions that involve a single atom are intramolecular"); + + /* Assign single-body interactions to the home zone. + * Position restraints are not handled here, but separately. + */ + if (iz == 0 && !(ftype == F_POSRES || ftype == F_FBPOSRES)) + { + bUse = true; + tiatoms[1] = atomIndex; + } + } + else if (nral == 2) + { + /* This is a two-body interaction, we can assign + * analogous to the non-bonded assignments. + */ + const int k_gl = (interactionConnectivity == InteractionConnectivity::Intramolecular) + ? + /* Get the global index using the offset in the molecule */ + (globalAtomIndex + iatoms[2] - atomIndexInMolecule) + : iatoms[2]; + if (const auto* entry = ga2la.find(k_gl)) + { + int kz = entry->cell; + if (kz >= zones->n) + { + kz -= zones->n; + } + /* Check zone interaction assignments */ + bUse = ((iz < iZones.ssize() && iz <= kz && iZones[iz].jZoneRange.isInRange(kz)) + || (kz < iZones.ssize() && iz > kz && iZones[kz].jZoneRange.isInRange(iz))); + if (bUse) + { + GMX_ASSERT(ftype != F_CONSTR || (iz == 0 && kz == 0), + "Constraint assigned here should only involve home atoms"); + + tiatoms[1] = atomIndex; + tiatoms[2] = entry->la; + /* If necessary check the cgcm distance */ + if (checkDistanceTwoBody + && dd_dist2(pbc_null, coordinates, tiatoms[1], tiatoms[2]) >= cutoffSquared) + { + bUse = false; + } + } + } + else + { + bUse = false; + } + } + else + { + /* Assign this multi-body bonded interaction to + * the local domain if we have all the atoms involved + * (local or communicated) and the minimum zone shift + * in each dimension is zero, for dimensions + * with 2 DD cells an extra check may be necessary. + */ + ivec k_zero, k_plus; + bUse = true; + clear_ivec(k_zero); + clear_ivec(k_plus); + for (int k = 1; k <= nral && bUse; k++) + { + const int k_gl = + (interactionConnectivity == InteractionConnectivity::Intramolecular) + ? + /* Get the global index using the offset in the molecule */ + (globalAtomIndex + iatoms[k] - atomIndexInMolecule) + : iatoms[k]; + const auto* entry = ga2la.find(k_gl); + if (entry == nullptr || entry->cell >= zones->n) + { + /* We do not have this atom of this interaction + * locally, or it comes from more than one cell + * away. + */ + bUse = FALSE; + } + else + { + tiatoms[k] = entry->la; + for (int d = 0; d < DIM; d++) + { + if (zones->shift[entry->cell][d] == 0) + { + k_zero[d] = k; + } + else + { + k_plus[d] = k; + } + } + } + } + bUse = (bUse && (k_zero[XX] != 0) && (k_zero[YY] != 0) && (k_zero[ZZ] != 0)); + if (checkDistanceMultiBody) + { + for (int d = 0; (d < DIM && bUse); d++) + { + /* Check if the distance falls within + * the cut-off to avoid possible multiple + * assignments of bonded interactions. + */ + if (rcheck[d] && k_plus[d] + && dd_dist2(pbc_null, coordinates, tiatoms[k_zero[d]], tiatoms[k_plus[d]]) + >= cutoffSquared) + { + bUse = FALSE; + } + } + } + } + if (bUse) + { + /* Add this interaction to the local topology */ + idef->il[ftype].push_back(tiatoms[0], nral, tiatoms + 1); + /* Sum so we can check in global_stat + * if we have everything. + */ + if (ddBondedChecking == DDBondedChecking::All + || !(interaction_function[ftype].flags & IF_LIMZERO)) + { + numBondedInteractions++; + } + } + } + j += 1 + nral_rt(ftype); + } + + return numBondedInteractions; +} + +/*! \brief Determine whether the local domain has responsibility for + * any of the position restraints for local atom \p atomIndex + * and assign those to the local domain. + * + * \returns The total number of bonded interactions for this atom for + * which this domain is responsible. + */ +static inline int assignPositionRestraintsForAtom(const int atomIndex, + const int mol, + const int atomIndexInMolecule, + const int numAtomsInMolecule, + gmx::ArrayRef rtil, + const int ind_start, + const int ind_end, + const gmx_molblock_t& molb, + const t_iparams* ip_in, + InteractionDefinitions* idef) +{ + constexpr int nral = 1; + + int numBondedInteractions = 0; + + int j = ind_start; + while (j < ind_end) + { + const int ftype = rtil[j++]; + auto iatoms = gmx::constArrayRefFromArray(rtil.data() + j, rtil.size() - j); + + if (ftype == F_POSRES || ftype == F_FBPOSRES) + { + std::array tiatoms = { iatoms[0], atomIndex }; + if (ftype == F_POSRES) + { + add_posres(mol, atomIndexInMolecule, numAtomsInMolecule, molb, tiatoms, ip_in, idef); + } + else + { + add_fbposres(mol, atomIndexInMolecule, numAtomsInMolecule, molb, tiatoms, ip_in, idef); + } + idef->il[ftype].push_back(tiatoms[0], nral, tiatoms.data() + 1); + numBondedInteractions++; + } + j += 1 + nral_rt(ftype); + } + + return numBondedInteractions; +} + +/*! \brief This function looks up and assigns bonded interactions for zone iz. + * + * With thread parallelizing each thread acts on a different atom range: + * at_start to at_end. + */ +static int make_bondeds_zone(gmx_reverse_top_t* rt, + ArrayRef globalAtomIndices, + const gmx_ga2la_t& ga2la, + const gmx_domdec_zones_t* zones, + const std::vector& molb, + const bool checkDistanceMultiBody, + const ivec rcheck, + const bool checkDistanceTwoBody, + const real cutoffSquared, + const t_pbc* pbc_null, + ArrayRef coordinates, + const t_iparams* ip_in, + InteractionDefinitions* idef, + int izone, + const gmx::Range& atomRange) +{ + int mb = 0; + int mt = 0; + int mol = 0; + int atomIndexInMolecule = 0; + + const auto ddBondedChecking = rt->options().ddBondedChecking; + + int numBondedInteractions = 0; + + for (int i : atomRange) + { + /* Get the global atom number */ + const int globalAtomIndex = globalAtomIndices[i]; + global_atomnr_to_moltype_ind( + rt->molblockIndices(), globalAtomIndex, &mb, &mt, &mol, &atomIndexInMolecule); + /* Check all intramolecular interactions assigned to this atom */ + gmx::ArrayRef index = rt->interactionListForMoleculeType(mt).index; + gmx::ArrayRef rtil = rt->interactionListForMoleculeType(mt).il; + + numBondedInteractions += assignInteractionsForAtom( + i, + globalAtomIndex, + atomIndexInMolecule, + index, + rtil, + index[atomIndexInMolecule], + index[atomIndexInMolecule + 1], + ga2la, + zones, + checkDistanceMultiBody, + rcheck, + checkDistanceTwoBody, + cutoffSquared, + pbc_null, + coordinates, + idef, + izone, + ddBondedChecking); + + // Assign position restraints, when present, for the home zone + if (izone == 0 && rt->hasPositionRestraints()) + { + numBondedInteractions += assignPositionRestraintsForAtom( + i, + mol, + atomIndexInMolecule, + rt->interactionListForMoleculeType(mt).numAtomsInMolecule, + rtil, + index[atomIndexInMolecule], + index[atomIndexInMolecule + 1], + molb[mb], + ip_in, + idef); + } + + if (rt->hasIntermolecularInteractions()) + { + /* Check all intermolecular interactions assigned to this atom */ + index = rt->interactionListForIntermolecularInteractions().index; + rtil = rt->interactionListForIntermolecularInteractions().il; + + numBondedInteractions += assignInteractionsForAtom( + i, + -1, // not used + -1, // not used + index, + rtil, + index[globalAtomIndex], + index[globalAtomIndex + 1], + ga2la, + zones, + checkDistanceMultiBody, + rcheck, + checkDistanceTwoBody, + cutoffSquared, + pbc_null, + coordinates, + idef, + izone, + ddBondedChecking); + } + } + + return numBondedInteractions; +} + +/*! \brief Set the exclusion data for i-zone \p iz */ +static void make_exclusions_zone(ArrayRef globalAtomIndices, + const gmx_ga2la_t& ga2la, + gmx_domdec_zones_t* zones, + ArrayRef molblockIndices, + const std::vector& moltype, + const int* atomInfo, + ListOfLists* lexcls, + int iz, + int at_start, + int at_end, + const gmx::ArrayRef intermolecularExclusionGroup) +{ + const auto& jAtomRange = zones->iZones[iz].jAtomRange; + + const gmx::index oldNumLists = lexcls->ssize(); + + std::vector exclusionsForAtom; + for (int at = at_start; at < at_end; at++) + { + exclusionsForAtom.clear(); + + if (GET_CGINFO_EXCL_INTER(atomInfo[at])) + { + int mb = 0; + int mt = 0; + int mol = 0; + int a_mol = 0; + + /* Copy the exclusions from the global top */ + int a_gl = globalAtomIndices[at]; + global_atomnr_to_moltype_ind(molblockIndices, a_gl, &mb, &mt, &mol, &a_mol); + const auto excls = moltype[mt].excls[a_mol]; + for (const int aj_mol : excls) + { + if (const auto* jEntry = ga2la.find(a_gl + aj_mol - a_mol)) + { + /* This check is not necessary, but it can reduce + * the number of exclusions in the list, which in turn + * can speed up the pair list construction a bit. + */ + if (jAtomRange.isInRange(jEntry->la)) + { + exclusionsForAtom.push_back(jEntry->la); + } + } + } + } + + bool isExcludedAtom = !intermolecularExclusionGroup.empty() + && std::find(intermolecularExclusionGroup.begin(), + intermolecularExclusionGroup.end(), + globalAtomIndices[at]) + != intermolecularExclusionGroup.end(); + + if (isExcludedAtom) + { + for (int qmAtomGlobalIndex : intermolecularExclusionGroup) + { + if (const auto* entry = ga2la.find(qmAtomGlobalIndex)) + { + exclusionsForAtom.push_back(entry->la); + } + } + } + + /* Append the exclusions for this atom to the topology */ + lexcls->pushBack(exclusionsForAtom); + } + + GMX_RELEASE_ASSERT( + lexcls->ssize() - oldNumLists == at_end - at_start, + "The number of exclusion list should match the number of atoms in the range"); +} + +/*! \brief Generate and store all required local bonded interactions in \p idef and local exclusions in \p lexcls + * + * \returns Total count of bonded interactions in the local topology on this domain */ +static int make_local_bondeds_excls(gmx_domdec_t* dd, + gmx_domdec_zones_t* zones, + const gmx_mtop_t& mtop, + const int* atomInfo, + const bool checkDistanceMultiBody, + const ivec rcheck, + const gmx_bool checkDistanceTwoBody, + const real cutoff, + const t_pbc* pbc_null, + ArrayRef coordinates, + InteractionDefinitions* idef, + ListOfLists* lexcls) +{ + int nzone_bondeds = 0; + + if (dd->reverse_top->hasInterAtomicInteractions()) + { + nzone_bondeds = zones->n; + } + else + { + /* Only single charge group (or atom) molecules, so interactions don't + * cross zone boundaries and we only need to assign in the home zone. + */ + nzone_bondeds = 1; + } + + /* We only use exclusions from i-zones to i- and j-zones */ + const int numIZonesForExclusions = (dd->haveExclusions ? zones->iZones.size() : 0); + + gmx_reverse_top_t* rt = dd->reverse_top.get(); + + const real cutoffSquared = gmx::square(cutoff); + + /* Clear the counts */ + idef->clear(); + int numBondedInteractions = 0; + + lexcls->clear(); + + for (int izone = 0; izone < nzone_bondeds; izone++) + { + const int cg0 = zones->cg_range[izone]; + const int cg1 = zones->cg_range[izone + 1]; + + gmx::ArrayRef threadWorkObjects = rt->threadWorkObjects(); + const int numThreads = threadWorkObjects.size(); +#pragma omp parallel for num_threads(numThreads) schedule(static) + for (int thread = 0; thread < numThreads; thread++) + { + try + { + InteractionDefinitions* idef_t = nullptr; + + int cg0t = cg0 + ((cg1 - cg0) * thread) / numThreads; + int cg1t = cg0 + ((cg1 - cg0) * (thread + 1)) / numThreads; + + if (thread == 0) + { + idef_t = idef; + } + else + { + idef_t = &threadWorkObjects[thread].idef; + idef_t->clear(); + } + + threadWorkObjects[thread].numBondedInteractions = + make_bondeds_zone(rt, + dd->globalAtomIndices, + *dd->ga2la, + zones, + mtop.molblock, + checkDistanceMultiBody, + rcheck, + checkDistanceTwoBody, + cutoffSquared, + pbc_null, + coordinates, + idef->iparams.data(), + idef_t, + izone, + gmx::Range(cg0t, cg1t)); + + if (izone < numIZonesForExclusions) + { + ListOfLists* excl_t = nullptr; + if (thread == 0) + { + // Thread 0 stores exclusions directly in the final storage + excl_t = lexcls; + } + else + { + // Threads > 0 store in temporary storage, starting at list index 0 + excl_t = &threadWorkObjects[thread].excl; + excl_t->clear(); + } + + /* No charge groups and no distance check required */ + make_exclusions_zone(dd->globalAtomIndices, + *dd->ga2la, + zones, + rt->molblockIndices(), + mtop.moltype, + atomInfo, + excl_t, + izone, + cg0t, + cg1t, + mtop.intermolecularExclusionGroup); + } + } + GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR + } + + if (threadWorkObjects.size() > 1) + { + combine_idef(idef, threadWorkObjects); + } + + for (const thread_work_t& th_work : threadWorkObjects) + { + numBondedInteractions += th_work.numBondedInteractions; + } + + if (izone < numIZonesForExclusions) + { + for (std::size_t th = 1; th < threadWorkObjects.size(); th++) + { + lexcls->appendListOfLists(threadWorkObjects[th].excl); + } + } + } + + if (debug) + { + fprintf(debug, "We have %d exclusions\n", lexcls->numElements()); + } + + return numBondedInteractions; +} + +int dd_make_local_top(gmx_domdec_t* dd, + gmx_domdec_zones_t* zones, + int npbcdim, + matrix box, + rvec cellsize_min, + const ivec npulse, + t_forcerec* fr, + ArrayRef coordinates, + const gmx_mtop_t& mtop, + gmx_localtop_t* ltop) +{ + real rc = -1; + ivec rcheck; + t_pbc pbc, *pbc_null = nullptr; + + if (debug) + { + fprintf(debug, "Making local topology\n"); + } + + bool checkDistanceMultiBody = false; + bool checkDistanceTwoBody = false; + + if (dd->reverse_top->hasInterAtomicInteractions()) + { + /* We need to check to which cell bondeds should be assigned */ + rc = dd_cutoff_twobody(dd); + if (debug) + { + fprintf(debug, "Two-body bonded cut-off distance is %g\n", rc); + } + + /* Should we check distances when assigning bonded interactions? */ + for (int d = 0; d < DIM; d++) + { + rcheck[d] = FALSE; + /* Only need to check for dimensions where the part of the box + * that is not communicated is smaller than the cut-off. + */ + if (d < npbcdim && dd->numCells[d] > 1 + && (dd->numCells[d] - npulse[d]) * cellsize_min[d] < 2 * rc) + { + if (dd->numCells[d] == 2) + { + rcheck[d] = TRUE; + checkDistanceMultiBody = TRUE; + } + /* Check for interactions between two atoms, + * where we can allow interactions up to the cut-off, + * instead of up to the smallest cell dimension. + */ + checkDistanceTwoBody = TRUE; + } + if (debug) + { + fprintf(debug, + "dim %d cellmin %f bonded rcheck[%d] = %d, checkDistanceTwoBody = %s\n", + d, + cellsize_min[d], + d, + rcheck[d], + gmx::boolToString(checkDistanceTwoBody)); + } + } + if (checkDistanceMultiBody || checkDistanceTwoBody) + { + if (fr->bMolPBC) + { + pbc_null = set_pbc_dd(&pbc, fr->pbcType, dd->numCells, TRUE, box); + } + else + { + pbc_null = nullptr; + } + } + } + + int numBondedInteractionsToReduce = make_local_bondeds_excls(dd, + zones, + mtop, + fr->atomInfo.data(), + checkDistanceMultiBody, + rcheck, + checkDistanceTwoBody, + rc, + pbc_null, + coordinates, + <op->idef, + <op->excls); + + /* The ilist is not sorted yet, + * we can only do this when we have the charge arrays. + */ + ltop->idef.ilsort = ilsortUNKNOWN; + + return numBondedInteractionsToReduce; +} + +void dd_sort_local_top(const gmx_domdec_t& dd, const t_mdatoms* mdatoms, gmx_localtop_t* ltop) +{ + if (dd.reverse_top->doSorting()) + { + gmx_sort_ilist_fe(<op->idef, mdatoms->chargeA, mdatoms->chargeB); + } + else + { + ltop->idef.ilsort = ilsortNO_FE; + } +} diff --git a/src/gromacs/domdec/localtopology.h b/src/gromacs/domdec/localtopology.h new file mode 100644 index 0000000000..6adfc26554 --- /dev/null +++ b/src/gromacs/domdec/localtopology.h @@ -0,0 +1,79 @@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright (c) 2021, 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. + * + * GROMACS is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + * + * GROMACS is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with GROMACS; if not, see + * http://www.gnu.org/licenses, or write to the Free Software Foundation, + * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + * + * If you want to redistribute modifications to GROMACS, please + * consider that scientific software is very special. Version + * control is crucial - bugs must be traceable. We will be happy to + * consider code for inclusion in the official distribution, but + * derived work must not be called official GROMACS. Details are found + * in the README & COPYING files - if they are missing, get the + * official version at http://www.gromacs.org. + * + * To help us fund GROMACS development, we humbly ask that you cite + * the research papers on the package. Check out http://www.gromacs.org. + */ +/*! \libinternal \file + * + * \brief This file makes declarations used for building + * the local topology + * + * \inlibraryapi + * \ingroup module_domdec + */ + +#ifndef GMX_DOMDEC_LOCALTOPOLOGY_H +#define GMX_DOMDEC_LOCALTOPOLOGY_H + +#include "gromacs/math/vectypes.h" + +struct gmx_domdec_t; +struct gmx_domdec_zones_t; +struct gmx_localtop_t; +struct gmx_mtop_t; +struct t_forcerec; +struct t_mdatoms; + +namespace gmx +{ +template +class ArrayRef; +} + +/*! \brief Generate the local topology and virtual site data + * + * \returns Total count of bonded interactions in the local topology on this domain */ +int dd_make_local_top(struct gmx_domdec_t* dd, + struct gmx_domdec_zones_t* zones, + int npbcdim, + matrix box, + rvec cellsize_min, + const ivec npulse, + t_forcerec* fr, + gmx::ArrayRef coordinates, + const gmx_mtop_t& top, + gmx_localtop_t* ltop); + +/*! \brief Sort ltop->ilist when we are doing free energy. */ +void dd_sort_local_top(const gmx_domdec_t& dd, const t_mdatoms* mdatoms, gmx_localtop_t* ltop); + +#endif diff --git a/src/gromacs/domdec/mdsetup.cpp b/src/gromacs/domdec/mdsetup.cpp index 47e4203f82..631d167227 100644 --- a/src/gromacs/domdec/mdsetup.cpp +++ b/src/gromacs/domdec/mdsetup.cpp @@ -38,6 +38,7 @@ #include "gromacs/domdec/domdec.h" #include "gromacs/domdec/domdec_struct.h" +#include "gromacs/domdec/localtopology.h" #include "gromacs/ewald/pme.h" #include "gromacs/listed_forces/listed_forces.h" #include "gromacs/mdlib/constr.h" diff --git a/src/gromacs/domdec/partition.cpp b/src/gromacs/domdec/partition.cpp index 2b1eea9948..dc826b55e7 100644 --- a/src/gromacs/domdec/partition.cpp +++ b/src/gromacs/domdec/partition.cpp @@ -56,10 +56,10 @@ #include "gromacs/domdec/collect.h" #include "gromacs/domdec/dlb.h" #include "gromacs/domdec/dlbtiming.h" -#include "gromacs/domdec/domdec.h" #include "gromacs/domdec/domdec_network.h" #include "gromacs/domdec/ga2la.h" #include "gromacs/domdec/localatomsetmanager.h" +#include "gromacs/domdec/localtopology.h" #include "gromacs/domdec/localtopologychecker.h" #include "gromacs/domdec/mdsetup.h" #include "gromacs/domdec/nsgrid.h" -- 2.22.0