src/gromacs/domdec/localtopology.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2006 - 2014, The GROMACS development team.
   5  * Copyright (c) 2015,2016,2017,2018,2019,2020,2021, by the GROMACS development team, led by
   6  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   7  * and including many others, as listed in the AUTHORS file in the
   8  * top-level source directory and at http://www.gromacs.org.
   9  *
  10  * GROMACS is free software; you can redistribute it and/or
  11  * modify it under the terms of the GNU Lesser General Public License
  12  * as published by the Free Software Foundation; either version 2.1
  13  * of the License, or (at your option) any later version.
  14  *
  15  * GROMACS is distributed in the hope that it will be useful,
  16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  18  * Lesser General Public License for more details.
  19  *
  20  * You should have received a copy of the GNU Lesser General Public
  21  * License along with GROMACS; if not, see
  22  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  23  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  24  *
  25  * If you want to redistribute modifications to GROMACS, please
  26  * consider that scientific software is very special. Version
  27  * control is crucial - bugs must be traceable. We will be happy to
  28  * consider code for inclusion in the official distribution, but
  29  * derived work must not be called official GROMACS. Details are found
  30  * in the README & COPYING files - if they are missing, get the
  31  * official version at http://www.gromacs.org.
  32  *
  33  * To help us fund GROMACS development, we humbly ask that you cite
  34  * the research papers on the package. Check out http://www.gromacs.org.
  35  */
  36
  37 /*! \internal \file
  38  *
  39  * \brief This file defines functions used by the domdec module while
  40  * building topologies local to a domain
  41  *
  42  * \author Berk Hess <hess@kth.se>
  43  * \ingroup module_domdec
  44  */
  45
  46 #include "gmxpre.h"
  47
  48 #include "gromacs/domdec/localtopology.h"
  49
  50 #include <algorithm>
  51 #include <iterator>
  52 #include <vector>
  53
  54 #include "gromacs/domdec/domdec_internal.h"
  55 #include "gromacs/domdec/domdec_struct.h"
  56 #include "gromacs/domdec/ga2la.h"
  57 #include "gromacs/domdec/options.h"
  58 #include "gromacs/domdec/reversetopology.h"
  59 #include "gromacs/math/vec.h"
  60 #include "gromacs/mdtypes/atominfo.h"
  61 #include "gromacs/mdtypes/forcerec.h"
  62 #include "gromacs/mdtypes/mdatom.h"
  63 #include "gromacs/pbcutil/pbc.h"
  64 #include "gromacs/topology/idef.h"
  65 #include "gromacs/topology/mtop_util.h"
  66 #include "gromacs/topology/topsort.h"
  67 #include "gromacs/utility/arrayref.h"
  68 #include "gromacs/utility/fatalerror.h"
  69 #include "gromacs/utility/gmxassert.h"
  70 #include "gromacs/utility/listoflists.h"
  71 #include "gromacs/utility/strconvert.h"
  72
  73 using gmx::ArrayRef;
  74 using gmx::DDBondedChecking;
  75 using gmx::ListOfLists;
  76 using gmx::RVec;
  77
  78 //! Container for returning molecule type and index information for an atom
  79 struct AtomIndexSet
  80 {
  81     int local;          //!< The local index
  82     int global;         //!< The global index
  83     int withinMolecule; //!< The index in the molecule this atom belongs to
  84 };
  85
  86 //! Container for returning molecule type and index information for an atom
  87 struct AtomInMolblock
  88 {
  89     int molblockIndex;       //!< The index into the molecule block array
  90     int moleculeType;        //!< The molecule type
  91     int moleculeIndex;       //!< The index of the molecule in the molecule block
  92     int atomIndexInMolecule; //!< The index of the atom in the molecule
  93 };
  94
  95 /*! \brief Return global topology molecule information for global atom index \p i_gl */
  96 static AtomInMolblock atomInMolblockFromGlobalAtomnr(ArrayRef<const MolblockIndices> molblockIndices,
  97                                                      const int globalAtomIndex)
  98 {
  99     // Find the molecule block whose range of global atom indices
 100     // includes globalAtomIndex, by being the first for which
 101     // globalAtomIndex is not greater than its end.
 102     auto molblockIt = std::partition_point(
 103             molblockIndices.begin(),
 104             molblockIndices.end(),
 105             [globalAtomIndex](const MolblockIndices& mbi) { return mbi.a_end <= globalAtomIndex; });
 106
 107     AtomInMolblock aim;
 108
 109     aim.molblockIndex = std::distance(molblockIndices.begin(), molblockIt);
 110     aim.moleculeType  = molblockIt->type;
 111     aim.moleculeIndex = (globalAtomIndex - molblockIt->a_start) / molblockIt->natoms_mol;
 112     aim.atomIndexInMolecule =
 113             (globalAtomIndex - molblockIt->a_start) - (aim.moleculeIndex) * molblockIt->natoms_mol;
 114
 115     return aim;
 116 }
 117
 118 /*! \brief Store a vsite at the end of \p il, returns an array with parameter and atom indices
 119  *
 120  * This routine is very similar to add_ifunc, but vsites interactions
 121  * have more work to do than other kinds of interactions, and the
 122  * complex way nral (and thus vector contents) depends on ftype
 123  * confuses static analysis tools unless we fuse the vsite
 124  * atom-indexing organization code with the ifunc-adding code, so that
 125  * they can see that nral is the same value.
 126  *
 127  * \param[in] ga2la  The global to local atom index mapping
 128  * \param[in] nral   The number of atoms involved in this vsite
 129  * \param[in] isLocalVsite  Whether all atoms involved in this vsite are local
 130  * \param[in] atomIndexSet  The atom index set for the virtual site
 131  * \param[in] iatoms     Indices, 0: interaction type, 1: vsite (unused), 2 ...: constructing atoms
 132  * \param[in,out] il     The interaction list to add this vsite to
 133  *
 134  * \returns an array with the parameter index and the NRAL atom indices
 135  */
 136 static inline ArrayRef<const int> add_ifunc_for_vsites(const gmx_ga2la_t&  ga2la,
 137                                                        const int           nral,
 138                                                        const bool          isLocalVsite,
 139                                                        const AtomIndexSet& atomIndexSet,
 140                                                        ArrayRef<const int> iatoms,
 141                                                        InteractionList*    il)
 142 {
 143     std::array<int, 1 + MAXATOMLIST> tiatoms;
 144
 145     /* Copy the type */
 146     tiatoms[0] = iatoms[0];
 147
 148     if (isLocalVsite)
 149     {
 150         /* We know the local index of the first atom */
 151         tiatoms[1] = atomIndexSet.local;
 152     }
 153     else
 154     {
 155         /* Convert later in make_local_vsites */
 156         tiatoms[1] = -atomIndexSet.global - 1;
 157     }
 158
 159     GMX_ASSERT(nral >= 2 && nral <= 5, "Invalid nral for vsites");
 160     for (int k = 2; k < 1 + nral; k++)
 161     {
 162         int ak_gl = atomIndexSet.global + iatoms[k] - atomIndexSet.withinMolecule;
 163         if (const int* homeIndex = ga2la.findHome(ak_gl))
 164         {
 165             tiatoms[k] = *homeIndex;
 166         }
 167         else
 168         {
 169             /* Copy the global index, convert later in make_local_vsites */
 170             tiatoms[k] = -(ak_gl + 1);
 171         }
 172         // Note that ga2la_get_home always sets the third parameter if
 173         // it returns TRUE
 174     }
 175     il->push_back(tiatoms[0], nral, tiatoms.data() + 1);
 176
 177     // Return an array with the parameter index and the nral atom indices
 178     return ArrayRef<int>(il->iatoms).subArray(il->iatoms.size() - (1 + nral), 1 + nral);
 179 }
 180
 181 /*! \brief Store a position restraint in idef and iatoms, complex because the parameters are different for each entry */
 182 static void add_posres(int                     mol,
 183                        int                     a_mol,
 184                        int                     numAtomsInMolecule,
 185                        const gmx_molblock_t&   molb,
 186                        gmx::ArrayRef<int>      iatoms,
 187                        const t_iparams*        ip_in,
 188                        InteractionDefinitions* idef)
 189 {
 190     /* This position restraint has not been added yet,
 191      * so it's index is the current number of position restraints.
 192      */
 193     const int n = idef->il[F_POSRES].size() / 2;
 194
 195     /* Get the position restraint coordinates from the molblock */
 196     const int a_molb = mol * numAtomsInMolecule + a_mol;
 197     GMX_ASSERT(a_molb < ssize(molb.posres_xA),
 198                "We need a sufficient number of position restraint coordinates");
 199     /* Copy the force constants */
 200     t_iparams ip        = ip_in[iatoms[0]];
 201     ip.posres.pos0A[XX] = molb.posres_xA[a_molb][XX];
 202     ip.posres.pos0A[YY] = molb.posres_xA[a_molb][YY];
 203     ip.posres.pos0A[ZZ] = molb.posres_xA[a_molb][ZZ];
 204     if (!molb.posres_xB.empty())
 205     {
 206         ip.posres.pos0B[XX] = molb.posres_xB[a_molb][XX];
 207         ip.posres.pos0B[YY] = molb.posres_xB[a_molb][YY];
 208         ip.posres.pos0B[ZZ] = molb.posres_xB[a_molb][ZZ];
 209     }
 210     else
 211     {
 212         ip.posres.pos0B[XX] = ip.posres.pos0A[XX];
 213         ip.posres.pos0B[YY] = ip.posres.pos0A[YY];
 214         ip.posres.pos0B[ZZ] = ip.posres.pos0A[ZZ];
 215     }
 216     /* Set the parameter index for idef->iparams_posres */
 217     iatoms[0] = n;
 218     idef->iparams_posres.push_back(ip);
 219
 220     GMX_ASSERT(int(idef->iparams_posres.size()) == n + 1,
 221                "The index of the parameter type should match n");
 222 }
 223
 224 /*! \brief Store a flat-bottomed position restraint in idef and iatoms, complex because the parameters are different for each entry */
 225 static void add_fbposres(int                     mol,
 226                          int                     a_mol,
 227                          int                     numAtomsInMolecule,
 228                          const gmx_molblock_t&   molb,
 229                          gmx::ArrayRef<int>      iatoms,
 230                          const t_iparams*        ip_in,
 231                          InteractionDefinitions* idef)
 232 {
 233     /* This flat-bottom position restraint has not been added yet,
 234      * so it's index is the current number of position restraints.
 235      */
 236     const int n = idef->il[F_FBPOSRES].size() / 2;
 237
 238     /* Get the position restraint coordinats from the molblock */
 239     const int a_molb = mol * numAtomsInMolecule + a_mol;
 240     GMX_ASSERT(a_molb < ssize(molb.posres_xA),
 241                "We need a sufficient number of position restraint coordinates");
 242     /* Copy the force constants */
 243     t_iparams ip = ip_in[iatoms[0]];
 244     /* Take reference positions from A position of normal posres */
 245     ip.fbposres.pos0[XX] = molb.posres_xA[a_molb][XX];
 246     ip.fbposres.pos0[YY] = molb.posres_xA[a_molb][YY];
 247     ip.fbposres.pos0[ZZ] = molb.posres_xA[a_molb][ZZ];
 248
 249     /* Note: no B-type for flat-bottom posres */
 250
 251     /* Set the parameter index for idef->iparams_fbposres */
 252     iatoms[0] = n;
 253     idef->iparams_fbposres.push_back(ip);
 254
 255     GMX_ASSERT(int(idef->iparams_fbposres.size()) == n + 1,
 256                "The index of the parameter type should match n");
 257 }
 258
 259 /*! \brief Store a virtual site interaction, complex because of PBC and recursion */
 260 static void add_vsite(const gmx_ga2la_t&      ga2la,
 261                       const reverse_ilist_t&  reverseIlist,
 262                       const int               ftype,
 263                       const int               nral,
 264                       const bool              isLocalVsite,
 265                       const AtomIndexSet&     atomIndexSet,
 266                       ArrayRef<const int>     iatoms,
 267                       InteractionDefinitions* idef)
 268 {
 269     /* Add this interaction to the local topology */
 270     ArrayRef<const int> tiatoms =
 271             add_ifunc_for_vsites(ga2la, nral, isLocalVsite, atomIndexSet, iatoms, &idef->il[ftype]);
 272
 273     if (iatoms[1 + nral])
 274     {
 275         /* Check for recursion */
 276         for (int k = 2; k < 1 + nral; k++)
 277         {
 278             if ((iatoms[1 + nral] & (2 << k)) && (tiatoms[k] < 0))
 279             {
 280                 /* This construction atoms is a vsite and not a home atom */
 281                 if (gmx_debug_at)
 282                 {
 283                     fprintf(debug,
 284                             "Constructing atom %d of vsite atom %d is a vsite and non-home\n",
 285                             iatoms[k] + 1,
 286                             atomIndexSet.withinMolecule + 1);
 287                 }
 288                 /* Find the vsite construction */
 289
 290                 /* Check all interactions assigned to this atom */
 291                 int j = reverseIlist.index[iatoms[k]];
 292                 while (j < reverseIlist.index[iatoms[k] + 1])
 293                 {
 294                     int ftype_r = reverseIlist.il[j++];
 295                     int nral_r  = NRAL(ftype_r);
 296                     if (interaction_function[ftype_r].flags & IF_VSITE)
 297                     {
 298                         /* Add this vsite (recursion) */
 299                         const AtomIndexSet atomIndexRecur = {
 300                             -1, atomIndexSet.global + iatoms[k] - iatoms[1], iatoms[k]
 301                         };
 302                         add_vsite(ga2la,
 303                                   reverseIlist,
 304                                   ftype_r,
 305                                   nral_r,
 306                                   false,
 307                                   atomIndexRecur,
 308                                   gmx::arrayRefFromArray(reverseIlist.il.data() + j,
 309                                                          reverseIlist.il.size() - j),
 310                                   idef);
 311                     }
 312                     j += 1 + nral_rt(ftype_r);
 313                 }
 314             }
 315         }
 316     }
 317 }
 318
 319 /*! \brief Returns the squared distance between atoms \p i and \p j */
 320 static real dd_dist2(const t_pbc* pbc_null, ArrayRef<const RVec> coordinates, const int i, const int j)
 321 {
 322     rvec dx;
 323
 324     if (pbc_null)
 325     {
 326         pbc_dx_aiuc(pbc_null, coordinates[i], coordinates[j], dx);
 327     }
 328     else
 329     {
 330         rvec_sub(coordinates[i], coordinates[j], dx);
 331     }
 332
 333     return norm2(dx);
 334 }
 335
 336 /*! \brief Append t_idef structures 1 to nsrc in src to *dest */
 337 static void combine_idef(InteractionDefinitions* dest, gmx::ArrayRef<const thread_work_t> src)
 338 {
 339     for (int ftype = 0; ftype < F_NRE; ftype++)
 340     {
 341         int n = 0;
 342         for (gmx::index s = 1; s < src.ssize(); s++)
 343         {
 344             n += src[s].idef.il[ftype].size();
 345         }
 346         if (n > 0)
 347         {
 348             for (gmx::index s = 1; s < src.ssize(); s++)
 349             {
 350                 dest->il[ftype].append(src[s].idef.il[ftype]);
 351             }
 352
 353             /* Position restraints need an additional treatment */
 354             if (ftype == F_POSRES || ftype == F_FBPOSRES)
 355             {
 356                 int                     nposres = dest->il[ftype].size() / 2;
 357                 std::vector<t_iparams>& iparams_dest =
 358                         (ftype == F_POSRES ? dest->iparams_posres : dest->iparams_fbposres);
 359
 360                 /* Set nposres to the number of original position restraints in dest */
 361                 for (gmx::index s = 1; s < src.ssize(); s++)
 362                 {
 363                     nposres -= src[s].idef.il[ftype].size() / 2;
 364                 }
 365
 366                 for (gmx::index s = 1; s < src.ssize(); s++)
 367                 {
 368                     const std::vector<t_iparams>& iparams_src =
 369                             (ftype == F_POSRES ? src[s].idef.iparams_posres : src[s].idef.iparams_fbposres);
 370                     iparams_dest.insert(iparams_dest.end(), iparams_src.begin(), iparams_src.end());
 371
 372                     /* Correct the indices into iparams_posres */
 373                     for (int i = 0; i < src[s].idef.il[ftype].size() / 2; i++)
 374                     {
 375                         /* Correct the index into iparams_posres */
 376                         dest->il[ftype].iatoms[nposres * 2] = nposres;
 377                         nposres++;
 378                     }
 379                 }
 380                 GMX_RELEASE_ASSERT(
 381                         int(iparams_dest.size()) == nposres,
 382                         "The number of parameters should match the number of restraints");
 383             }
 384         }
 385     }
 386 }
 387
 388 /*! \brief Determine whether the local domain has responsibility for
 389  * any of the bonded interactions for local atom \p atomIndex
 390  * and assign those to the local domain.
 391  *
 392  * \returns The total number of bonded interactions for this atom for
 393  * which this domain is responsible.
 394  */
 395 static inline int assignInteractionsForAtom(const AtomIndexSet&       atomIndexSet,
 396                                             const reverse_ilist_t&    reverseIlist,
 397                                             const gmx_ga2la_t&        ga2la,
 398                                             const gmx_domdec_zones_t& zones,
 399                                             const bool                checkDistanceMultiBody,
 400                                             const ivec                rcheck,
 401                                             const bool                checkDistanceTwoBody,
 402                                             const real                cutoffSquared,
 403                                             const t_pbc*              pbc_null,
 404                                             ArrayRef<const RVec>      coordinates,
 405                                             InteractionDefinitions*   idef,
 406                                             const int                 iz,
 407                                             const DDBondedChecking    ddBondedChecking)
 408 {
 409     gmx::ArrayRef<const DDPairInteractionRanges> iZones = zones.iZones;
 410
 411     ArrayRef<const int> rtil = reverseIlist.il;
 412
 413     int numBondedInteractions = 0;
 414
 415     int       j        = reverseIlist.index[atomIndexSet.withinMolecule];
 416     const int indexEnd = reverseIlist.index[atomIndexSet.withinMolecule + 1];
 417     while (j < indexEnd)
 418     {
 419         int tiatoms[1 + MAXATOMLIST];
 420
 421         const int ftype  = rtil[j++];
 422         auto      iatoms = gmx::constArrayRefFromArray(rtil.data() + j, rtil.size() - j);
 423         const int nral   = NRAL(ftype);
 424         if (interaction_function[ftype].flags & IF_VSITE)
 425         {
 426             /* The vsite construction goes where the vsite itself is */
 427             if (iz == 0)
 428             {
 429                 add_vsite(ga2la, reverseIlist, ftype, nral, true, atomIndexSet, iatoms, idef);
 430             }
 431         }
 432         else
 433         {
 434             bool bUse = false;
 435
 436             /* Copy the type */
 437             tiatoms[0] = iatoms[0];
 438
 439             if (nral == 1)
 440             {
 441                 /* Assign single-body interactions to the home zone.
 442                  * Position restraints are not handled here, but separately.
 443                  */
 444                 if (iz == 0 && !(ftype == F_POSRES || ftype == F_FBPOSRES))
 445                 {
 446                     bUse       = true;
 447                     tiatoms[1] = atomIndexSet.local;
 448                 }
 449             }
 450             else if (nral == 2)
 451             {
 452                 /* This is a two-body interaction, we can assign
 453                  * analogous to the non-bonded assignments.
 454                  */
 455                 /* Get the global index using the offset in the molecule */
 456                 const int k_gl = atomIndexSet.global + iatoms[2] - atomIndexSet.withinMolecule;
 457                 if (const auto* entry = ga2la.find(k_gl))
 458                 {
 459                     int kz = entry->cell;
 460                     if (kz >= zones.n)
 461                     {
 462                         kz -= zones.n;
 463                     }
 464                     /* Check zone interaction assignments */
 465                     bUse = ((iz < iZones.ssize() && iz <= kz && iZones[iz].jZoneRange.isInRange(kz))
 466                             || (kz < iZones.ssize() && iz > kz && iZones[kz].jZoneRange.isInRange(iz)));
 467                     if (bUse)
 468                     {
 469                         GMX_ASSERT(ftype != F_CONSTR || (iz == 0 && kz == 0),
 470                                    "Constraint assigned here should only involve home atoms");
 471
 472                         tiatoms[1] = atomIndexSet.local;
 473                         tiatoms[2] = entry->la;
 474                         /* If necessary check the cgcm distance */
 475                         if (checkDistanceTwoBody
 476                             && dd_dist2(pbc_null, coordinates, tiatoms[1], tiatoms[2]) >= cutoffSquared)
 477                         {
 478                             bUse = false;
 479                         }
 480                     }
 481                 }
 482                 else
 483                 {
 484                     bUse = false;
 485                 }
 486             }
 487             else
 488             {
 489                 /* Assign this multi-body bonded interaction to
 490                  * the local domain if we have all the atoms involved
 491                  * (local or communicated) and the minimum zone shift
 492                  * in each dimension is zero, for dimensions
 493                  * with 2 DD cells an extra check may be necessary.
 494                  */
 495                 ivec k_zero, k_plus;
 496                 bUse = true;
 497                 clear_ivec(k_zero);
 498                 clear_ivec(k_plus);
 499                 for (int k = 1; k <= nral && bUse; k++)
 500                 {
 501                     /* Get the global index using the offset in the molecule */
 502                     const int k_gl = atomIndexSet.global + iatoms[k] - atomIndexSet.withinMolecule;
 503                     const auto* entry = ga2la.find(k_gl);
 504                     if (entry == nullptr || entry->cell >= zones.n)
 505                     {
 506                         /* We do not have this atom of this interaction
 507                          * locally, or it comes from more than one cell
 508                          * away.
 509                          */
 510                         bUse = FALSE;
 511                     }
 512                     else
 513                     {
 514                         tiatoms[k] = entry->la;
 515                         for (int d = 0; d < DIM; d++)
 516                         {
 517                             if (zones.shift[entry->cell][d] == 0)
 518                             {
 519                                 k_zero[d] = k;
 520                             }
 521                             else
 522                             {
 523                                 k_plus[d] = k;
 524                             }
 525                         }
 526                     }
 527                 }
 528                 bUse = (bUse && (k_zero[XX] != 0) && (k_zero[YY] != 0) && (k_zero[ZZ] != 0));
 529                 if (checkDistanceMultiBody)
 530                 {
 531                     for (int d = 0; (d < DIM && bUse); d++)
 532                     {
 533                         /* Check if the distance falls within
 534                          * the cut-off to avoid possible multiple
 535                          * assignments of bonded interactions.
 536                          */
 537                         if (rcheck[d] && k_plus[d]
 538                             && dd_dist2(pbc_null, coordinates, tiatoms[k_zero[d]], tiatoms[k_plus[d]])
 539                                        >= cutoffSquared)
 540                         {
 541                             bUse = FALSE;
 542                         }
 543                     }
 544                 }
 545             }
 546             if (bUse)
 547             {
 548                 /* Add this interaction to the local topology */
 549                 idef->il[ftype].push_back(tiatoms[0], nral, tiatoms + 1);
 550                 /* Sum so we can check in global_stat
 551                  * if we have everything.
 552                  */
 553                 if (ddBondedChecking == DDBondedChecking::All
 554                     || !(interaction_function[ftype].flags & IF_LIMZERO))
 555                 {
 556                     numBondedInteractions++;
 557                 }
 558             }
 559         }
 560         j += 1 + nral_rt(ftype);
 561     }
 562
 563     return numBondedInteractions;
 564 }
 565
 566 /*! \brief Determine whether the local domain has responsibility for
 567  * any of the position restraints for local atom \p atomIndex
 568  * and assign those to the local domain.
 569  *
 570  * \returns The total number of bonded interactions for this atom for
 571  * which this domain is responsible.
 572  */
 573 static inline int assignPositionRestraintsForAtom(const AtomIndexSet&     atomIndexSet,
 574                                                   const int               moleculeIndex,
 575                                                   const int               numAtomsInMolecule,
 576                                                   const reverse_ilist_t&  reverseIlist,
 577                                                   const gmx_molblock_t&   molb,
 578                                                   const t_iparams*        ip_in,
 579                                                   InteractionDefinitions* idef)
 580 {
 581     constexpr int nral = 1;
 582
 583     ArrayRef<const int> rtil = reverseIlist.il;
 584
 585     int numBondedInteractions = 0;
 586
 587     int       j        = reverseIlist.index[atomIndexSet.withinMolecule];
 588     const int indexEnd = reverseIlist.index[atomIndexSet.withinMolecule + 1];
 589     while (j < indexEnd)
 590     {
 591         const int ftype  = rtil[j++];
 592         auto      iatoms = gmx::constArrayRefFromArray(rtil.data() + j, rtil.size() - j);
 593
 594         if (ftype == F_POSRES || ftype == F_FBPOSRES)
 595         {
 596             std::array<int, 1 + nral> tiatoms = { iatoms[0], atomIndexSet.local };
 597             if (ftype == F_POSRES)
 598             {
 599                 add_posres(moleculeIndex, atomIndexSet.withinMolecule, numAtomsInMolecule, molb, tiatoms, ip_in, idef);
 600             }
 601             else
 602             {
 603                 add_fbposres(
 604                         moleculeIndex, atomIndexSet.withinMolecule, numAtomsInMolecule, molb, tiatoms, ip_in, idef);
 605             }
 606             idef->il[ftype].push_back(tiatoms[0], nral, tiatoms.data() + 1);
 607             numBondedInteractions++;
 608         }
 609         j += 1 + nral_rt(ftype);
 610     }
 611
 612     return numBondedInteractions;
 613 }
 614
 615 /*! \brief This function looks up and assigns bonded interactions for zone iz.
 616  *
 617  * With thread parallelizing each thread acts on a different atom range:
 618  * at_start to at_end.
 619  */
 620 static int make_bondeds_zone(const gmx_reverse_top_t&           rt,
 621                              ArrayRef<const int>                globalAtomIndices,
 622                              const gmx_ga2la_t&                 ga2la,
 623                              const gmx_domdec_zones_t&          zones,
 624                              const std::vector<gmx_molblock_t>& molb,
 625                              const bool                         checkDistanceMultiBody,
 626                              const ivec                         rcheck,
 627                              const bool                         checkDistanceTwoBody,
 628                              const real                         cutoffSquared,
 629                              const t_pbc*                       pbc_null,
 630                              ArrayRef<const RVec>               coordinates,
 631                              const t_iparams*                   ip_in,
 632                              InteractionDefinitions*            idef,
 633                              int                                izone,
 634                              const gmx::Range<int>&             atomRange)
 635 {
 636     const auto ddBondedChecking = rt.options().ddBondedChecking;
 637
 638     int numBondedInteractions = 0;
 639
 640     for (int atomIndexLocal : atomRange)
 641     {
 642         /* Get the global atom number */
 643         const int  atomIndexGlobal = globalAtomIndices[atomIndexLocal];
 644         const auto aim = atomInMolblockFromGlobalAtomnr(rt.molblockIndices(), atomIndexGlobal);
 645
 646         const AtomIndexSet atomIndexMol = { atomIndexLocal, atomIndexGlobal, aim.atomIndexInMolecule };
 647         const auto&        ilistMol     = rt.interactionListForMoleculeType(aim.moleculeType);
 648         numBondedInteractions += assignInteractionsForAtom(atomIndexMol,
 649                                                            ilistMol,
 650                                                            ga2la,
 651                                                            zones,
 652                                                            checkDistanceMultiBody,
 653                                                            rcheck,
 654                                                            checkDistanceTwoBody,
 655                                                            cutoffSquared,
 656                                                            pbc_null,
 657                                                            coordinates,
 658                                                            idef,
 659                                                            izone,
 660                                                            ddBondedChecking);
 661
 662         // Assign position restraints, when present, for the home zone
 663         if (izone == 0 && rt.hasPositionRestraints())
 664         {
 665             numBondedInteractions +=
 666                     assignPositionRestraintsForAtom(atomIndexMol,
 667                                                     aim.moleculeIndex,
 668                                                     ilistMol.numAtomsInMolecule,
 669                                                     rt.interactionListForMoleculeType(aim.moleculeType),
 670                                                     molb[aim.molblockIndex],
 671                                                     ip_in,
 672                                                     idef);
 673         }
 674
 675         if (rt.hasIntermolecularInteractions())
 676         {
 677             /* Check all intermolecular interactions assigned to this atom.
 678              * Note that we will index the intermolecular reverse ilist with atomIndexGlobal.
 679              */
 680             const AtomIndexSet atomIndexIntermol = { atomIndexLocal, atomIndexGlobal, atomIndexGlobal };
 681             numBondedInteractions +=
 682                     assignInteractionsForAtom(atomIndexIntermol,
 683                                               rt.interactionListForIntermolecularInteractions(),
 684                                               ga2la,
 685                                               zones,
 686                                               checkDistanceMultiBody,
 687                                               rcheck,
 688                                               checkDistanceTwoBody,
 689                                               cutoffSquared,
 690                                               pbc_null,
 691                                               coordinates,
 692                                               idef,
 693                                               izone,
 694                                               ddBondedChecking);
 695         }
 696     }
 697
 698     return numBondedInteractions;
 699 }
 700
 701 /*! \brief Set the exclusion data for i-zone \p iz */
 702 static void make_exclusions_zone(ArrayRef<const int>               globalAtomIndices,
 703                                  const gmx_ga2la_t&                ga2la,
 704                                  const gmx_domdec_zones_t&         zones,
 705                                  ArrayRef<const MolblockIndices>   molblockIndices,
 706                                  const std::vector<gmx_moltype_t>& moltype,
 707                                  gmx::ArrayRef<const int64_t>      atomInfo,
 708                                  ListOfLists<int>*                 lexcls,
 709                                  int                               iz,
 710                                  int                               at_start,
 711                                  int                               at_end,
 712                                  const gmx::ArrayRef<const int>    intermolecularExclusionGroup)
 713 {
 714     const auto& jAtomRange = zones.iZones[iz].jAtomRange;
 715
 716     const gmx::index oldNumLists = lexcls->ssize();
 717
 718     std::vector<int> exclusionsForAtom;
 719     for (int at = at_start; at < at_end; at++)
 720     {
 721         exclusionsForAtom.clear();
 722
 723         if (atomInfo[at] & gmx::sc_atomInfo_Exclusion)
 724         {
 725             /* Copy the exclusions from the global top */
 726             const int                          globalAtomIndex = globalAtomIndices[at];
 727             const MolecularTopologyAtomIndices mtai =
 728                     globalAtomIndexToMoltypeIndices(molblockIndices, globalAtomIndex);
 729             const auto& excls = moltype[mtai.moleculeType].excls[mtai.atomIndex];
 730             for (const int excludedAtomIndexInMolecule : excls)
 731             {
 732                 if (const auto* jEntry =
 733                             ga2la.find(globalAtomIndex + excludedAtomIndexInMolecule - mtai.atomIndex))
 734                 {
 735                     /* This check is not necessary, but it can reduce
 736                      * the number of exclusions in the list, which in turn
 737                      * can speed up the pair list construction a bit.
 738                      */
 739                     if (jAtomRange.isInRange(jEntry->la))
 740                     {
 741                         exclusionsForAtom.push_back(jEntry->la);
 742                     }
 743                 }
 744             }
 745         }
 746
 747         bool isExcludedAtom = !intermolecularExclusionGroup.empty()
 748                               && std::find(intermolecularExclusionGroup.begin(),
 749                                            intermolecularExclusionGroup.end(),
 750                                            globalAtomIndices[at])
 751                                          != intermolecularExclusionGroup.end();
 752
 753         if (isExcludedAtom)
 754         {
 755             for (int qmAtomGlobalIndex : intermolecularExclusionGroup)
 756             {
 757                 if (const auto* entry = ga2la.find(qmAtomGlobalIndex))
 758                 {
 759                     exclusionsForAtom.push_back(entry->la);
 760                 }
 761             }
 762         }
 763
 764         /* Append the exclusions for this atom to the topology */
 765         lexcls->pushBack(exclusionsForAtom);
 766     }
 767
 768     GMX_RELEASE_ASSERT(
 769             lexcls->ssize() - oldNumLists == at_end - at_start,
 770             "The number of exclusion list should match the number of atoms in the range");
 771 }
 772
 773 /*! \brief Generate and store all required local bonded interactions in \p idef and local exclusions in \p lexcls
 774  *
 775  * \returns Total count of bonded interactions in the local topology on this domain */
 776 static int make_local_bondeds_excls(const gmx_domdec_t&       dd,
 777                                     const gmx_domdec_zones_t& zones,
 778                                     const gmx_mtop_t&         mtop,
 779                                     ArrayRef<const int64_t>   atomInfo,
 780                                     const bool                checkDistanceMultiBody,
 781                                     const ivec                rcheck,
 782                                     const gmx_bool            checkDistanceTwoBody,
 783                                     const real                cutoff,
 784                                     const t_pbc*              pbc_null,
 785                                     ArrayRef<const RVec>      coordinates,
 786                                     InteractionDefinitions*   idef,
 787                                     ListOfLists<int>*         lexcls)
 788 {
 789     int nzone_bondeds = 0;
 790
 791     if (dd.reverse_top->hasInterAtomicInteractions())
 792     {
 793         nzone_bondeds = zones.n;
 794     }
 795     else
 796     {
 797         /* Only single charge group (or atom) molecules, so interactions don't
 798          * cross zone boundaries and we only need to assign in the home zone.
 799          */
 800         nzone_bondeds = 1;
 801     }
 802
 803     /* We only use exclusions from i-zones to i- and j-zones */
 804     const int numIZonesForExclusions = (dd.haveExclusions ? zones.iZones.size() : 0);
 805
 806     const gmx_reverse_top_t& rt = *dd.reverse_top;
 807
 808     const real cutoffSquared = gmx::square(cutoff);
 809
 810     /* Clear the counts */
 811     idef->clear();
 812     int numBondedInteractions = 0;
 813
 814     lexcls->clear();
 815
 816     for (int izone = 0; izone < nzone_bondeds; izone++)
 817     {
 818         const int cg0 = zones.cg_range[izone];
 819         const int cg1 = zones.cg_range[izone + 1];
 820
 821         gmx::ArrayRef<thread_work_t> threadWorkObjects = rt.threadWorkObjects();
 822         const int                    numThreads        = threadWorkObjects.size();
 823 #pragma omp parallel for num_threads(numThreads) schedule(static)
 824         for (int thread = 0; thread < numThreads; thread++)
 825         {
 826             try
 827             {
 828                 InteractionDefinitions* idef_t = nullptr;
 829
 830                 int cg0t = cg0 + ((cg1 - cg0) * thread) / numThreads;
 831                 int cg1t = cg0 + ((cg1 - cg0) * (thread + 1)) / numThreads;
 832
 833                 if (thread == 0)
 834                 {
 835                     idef_t = idef;
 836                 }
 837                 else
 838                 {
 839                     idef_t = &threadWorkObjects[thread].idef;
 840                     idef_t->clear();
 841                 }
 842
 843                 threadWorkObjects[thread].numBondedInteractions =
 844                         make_bondeds_zone(rt,
 845                                           dd.globalAtomIndices,
 846                                           *dd.ga2la,
 847                                           zones,
 848                                           mtop.molblock,
 849                                           checkDistanceMultiBody,
 850                                           rcheck,
 851                                           checkDistanceTwoBody,
 852                                           cutoffSquared,
 853                                           pbc_null,
 854                                           coordinates,
 855                                           idef->iparams.data(),
 856                                           idef_t,
 857                                           izone,
 858                                           gmx::Range<int>(cg0t, cg1t));
 859
 860                 if (izone < numIZonesForExclusions)
 861                 {
 862                     ListOfLists<int>* excl_t = nullptr;
 863                     if (thread == 0)
 864                     {
 865                         // Thread 0 stores exclusions directly in the final storage
 866                         excl_t = lexcls;
 867                     }
 868                     else
 869                     {
 870                         // Threads > 0 store in temporary storage, starting at list index 0
 871                         excl_t = &threadWorkObjects[thread].excl;
 872                         excl_t->clear();
 873                     }
 874
 875                     /* No charge groups and no distance check required */
 876                     make_exclusions_zone(dd.globalAtomIndices,
 877                                          *dd.ga2la,
 878                                          zones,
 879                                          rt.molblockIndices(),
 880                                          mtop.moltype,
 881                                          atomInfo,
 882                                          excl_t,
 883                                          izone,
 884                                          cg0t,
 885                                          cg1t,
 886                                          mtop.intermolecularExclusionGroup);
 887                 }
 888             }
 889             GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
 890         }
 891
 892         if (threadWorkObjects.size() > 1)
 893         {
 894             combine_idef(idef, threadWorkObjects);
 895         }
 896
 897         for (const thread_work_t& th_work : threadWorkObjects)
 898         {
 899             numBondedInteractions += th_work.numBondedInteractions;
 900         }
 901
 902         if (izone < numIZonesForExclusions)
 903         {
 904             for (std::size_t th = 1; th < threadWorkObjects.size(); th++)
 905             {
 906                 lexcls->appendListOfLists(threadWorkObjects[th].excl);
 907             }
 908         }
 909     }
 910
 911     if (debug)
 912     {
 913         fprintf(debug, "We have %d exclusions\n", lexcls->numElements());
 914     }
 915
 916     return numBondedInteractions;
 917 }
 918
 919 int dd_make_local_top(const gmx_domdec_t&          dd,
 920                       const gmx_domdec_zones_t&    zones,
 921                       int                          npbcdim,
 922                       matrix                       box,
 923                       rvec                         cellsize_min,
 924                       const ivec                   npulse,
 925                       t_forcerec*                  fr,
 926                       ArrayRef<const RVec>         coordinates,
 927                       const gmx_mtop_t&            mtop,
 928                       gmx::ArrayRef<const int64_t> atomInfo,
 929                       gmx_localtop_t*              ltop)
 930 {
 931     real  rc = -1;
 932     ivec  rcheck;
 933     t_pbc pbc, *pbc_null = nullptr;
 934
 935     if (debug)
 936     {
 937         fprintf(debug, "Making local topology\n");
 938     }
 939
 940     bool checkDistanceMultiBody = false;
 941     bool checkDistanceTwoBody   = false;
 942
 943     if (dd.reverse_top->hasInterAtomicInteractions())
 944     {
 945         /* We need to check to which cell bondeds should be assigned */
 946         rc = dd_cutoff_twobody(&dd);
 947         if (debug)
 948         {
 949             fprintf(debug, "Two-body bonded cut-off distance is %g\n", rc);
 950         }
 951
 952         /* Should we check distances when assigning bonded interactions? */
 953         for (int d = 0; d < DIM; d++)
 954         {
 955             rcheck[d] = FALSE;
 956             /* Only need to check for dimensions where the part of the box
 957              * that is not communicated is smaller than the cut-off.
 958              */
 959             if (d < npbcdim && dd.numCells[d] > 1 && (dd.numCells[d] - npulse[d]) * cellsize_min[d] < 2 * rc)
 960             {
 961                 if (dd.numCells[d] == 2)
 962                 {
 963                     rcheck[d]              = TRUE;
 964                     checkDistanceMultiBody = TRUE;
 965                 }
 966                 /* Check for interactions between two atoms,
 967                  * where we can allow interactions up to the cut-off,
 968                  * instead of up to the smallest cell dimension.
 969                  */
 970                 checkDistanceTwoBody = TRUE;
 971             }
 972             if (debug)
 973             {
 974                 fprintf(debug,
 975                         "dim %d cellmin %f bonded rcheck[%d] = %d, checkDistanceTwoBody = %s\n",
 976                         d,
 977                         cellsize_min[d],
 978                         d,
 979                         rcheck[d],
 980                         gmx::boolToString(checkDistanceTwoBody));
 981             }
 982         }
 983         if (checkDistanceMultiBody || checkDistanceTwoBody)
 984         {
 985             if (fr->bMolPBC)
 986             {
 987                 pbc_null = set_pbc_dd(&pbc, fr->pbcType, dd.numCells, TRUE, box);
 988             }
 989             else
 990             {
 991                 pbc_null = nullptr;
 992             }
 993         }
 994     }
 995
 996     int numBondedInteractionsToReduce = make_local_bondeds_excls(dd,
 997                                                                  zones,
 998                                                                  mtop,
 999                                                                  fr->atomInfo,
1000                                                                  checkDistanceMultiBody,
1001                                                                  rcheck,
1002                                                                  checkDistanceTwoBody,
1003                                                                  rc,
1004                                                                  pbc_null,
1005                                                                  coordinates,
1006                                                                  &ltop->idef,
1007                                                                  &ltop->excls);
1008
1009     if (dd.reverse_top->doListedForcesSorting())
1010     {
1011         gmx_sort_ilist_fe(&ltop->idef, atomInfo);
1012     }
1013     else
1014     {
1015         ltop->idef.ilsort = ilsortNO_FE;
1016     }
1017
1018     return numBondedInteractionsToReduce;
1019 }