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