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     int mb                  = 0;
 618     int mt                  = 0;
 619     int mol                 = 0;
 620     int atomIndexInMolecule = 0;
 621
 622     const auto ddBondedChecking = rt->options().ddBondedChecking;
 623
 624     int numBondedInteractions = 0;
 625
 626     for (int i : atomRange)
 627     {
 628         /* Get the global atom number */
 629         const int globalAtomIndex = globalAtomIndices[i];
 630         global_atomnr_to_moltype_ind(
 631                 rt->molblockIndices(), globalAtomIndex, &mb, &mt, &mol, &atomIndexInMolecule);
 632         /* Check all intramolecular interactions assigned to this atom */
 633         gmx::ArrayRef<const int>     index = rt->interactionListForMoleculeType(mt).index;
 634         gmx::ArrayRef<const t_iatom> rtil  = rt->interactionListForMoleculeType(mt).il;
 635
 636         numBondedInteractions += assignInteractionsForAtom<InteractionConnectivity::Intramolecular>(
 637                 i,
 638                 globalAtomIndex,
 639                 atomIndexInMolecule,
 640                 index,
 641                 rtil,
 642                 index[atomIndexInMolecule],
 643                 index[atomIndexInMolecule + 1],
 644                 ga2la,
 645                 zones,
 646                 checkDistanceMultiBody,
 647                 rcheck,
 648                 checkDistanceTwoBody,
 649                 cutoffSquared,
 650                 pbc_null,
 651                 coordinates,
 652                 idef,
 653                 izone,
 654                 ddBondedChecking);
 655
 656         // Assign position restraints, when present, for the home zone
 657         if (izone == 0 && rt->hasPositionRestraints())
 658         {
 659             numBondedInteractions += assignPositionRestraintsForAtom(
 660                     i,
 661                     mol,
 662                     atomIndexInMolecule,
 663                     rt->interactionListForMoleculeType(mt).numAtomsInMolecule,
 664                     rtil,
 665                     index[atomIndexInMolecule],
 666                     index[atomIndexInMolecule + 1],
 667                     molb[mb],
 668                     ip_in,
 669                     idef);
 670         }
 671
 672         if (rt->hasIntermolecularInteractions())
 673         {
 674             /* Check all intermolecular interactions assigned to this atom */
 675             index = rt->interactionListForIntermolecularInteractions().index;
 676             rtil  = rt->interactionListForIntermolecularInteractions().il;
 677
 678             numBondedInteractions += assignInteractionsForAtom<InteractionConnectivity::Intermolecular>(
 679                     i,
 680                     -1, // not used
 681                     -1, // not used
 682                     index,
 683                     rtil,
 684                     index[globalAtomIndex],
 685                     index[globalAtomIndex + 1],
 686                     ga2la,
 687                     zones,
 688                     checkDistanceMultiBody,
 689                     rcheck,
 690                     checkDistanceTwoBody,
 691                     cutoffSquared,
 692                     pbc_null,
 693                     coordinates,
 694                     idef,
 695                     izone,
 696                     ddBondedChecking);
 697         }
 698     }
 699
 700     return numBondedInteractions;
 701 }
 702
 703 /*! \brief Set the exclusion data for i-zone \p iz */
 704 static void make_exclusions_zone(ArrayRef<const int>               globalAtomIndices,
 705                                  const gmx_ga2la_t&                ga2la,
 706                                  gmx_domdec_zones_t*               zones,
 707                                  ArrayRef<const MolblockIndices>   molblockIndices,
 708                                  const std::vector<gmx_moltype_t>& moltype,
 709                                  gmx::ArrayRef<const int64_t>      atomInfo,
 710                                  ListOfLists<int>*                 lexcls,
 711                                  int                               iz,
 712                                  int                               at_start,
 713                                  int                               at_end,
 714                                  const gmx::ArrayRef<const int>    intermolecularExclusionGroup)
 715 {
 716     const auto& jAtomRange = zones->iZones[iz].jAtomRange;
 717
 718     const gmx::index oldNumLists = lexcls->ssize();
 719
 720     std::vector<int> exclusionsForAtom;
 721     for (int at = at_start; at < at_end; at++)
 722     {
 723         exclusionsForAtom.clear();
 724
 725         if (atomInfo[at] & gmx::sc_atomInfo_Exclusion)
 726         {
 727             int mb    = 0;
 728             int mt    = 0;
 729             int mol   = 0;
 730             int a_mol = 0;
 731
 732             /* Copy the exclusions from the global top */
 733             int a_gl = globalAtomIndices[at];
 734             global_atomnr_to_moltype_ind(molblockIndices, a_gl, &mb, &mt, &mol, &a_mol);
 735             const auto excls = moltype[mt].excls[a_mol];
 736             for (const int aj_mol : excls)
 737             {
 738                 if (const auto* jEntry = ga2la.find(a_gl + aj_mol - a_mol))
 739                 {
 740                     /* This check is not necessary, but it can reduce
 741                      * the number of exclusions in the list, which in turn
 742                      * can speed up the pair list construction a bit.
 743                      */
 744                     if (jAtomRange.isInRange(jEntry->la))
 745                     {
 746                         exclusionsForAtom.push_back(jEntry->la);
 747                     }
 748                 }
 749             }
 750         }
 751
 752         bool isExcludedAtom = !intermolecularExclusionGroup.empty()
 753                               && std::find(intermolecularExclusionGroup.begin(),
 754                                            intermolecularExclusionGroup.end(),
 755                                            globalAtomIndices[at])
 756                                          != intermolecularExclusionGroup.end();
 757
 758         if (isExcludedAtom)
 759         {
 760             for (int qmAtomGlobalIndex : intermolecularExclusionGroup)
 761             {
 762                 if (const auto* entry = ga2la.find(qmAtomGlobalIndex))
 763                 {
 764                     exclusionsForAtom.push_back(entry->la);
 765                 }
 766             }
 767         }
 768
 769         /* Append the exclusions for this atom to the topology */
 770         lexcls->pushBack(exclusionsForAtom);
 771     }
 772
 773     GMX_RELEASE_ASSERT(
 774             lexcls->ssize() - oldNumLists == at_end - at_start,
 775             "The number of exclusion list should match the number of atoms in the range");
 776 }
 777
 778 /*! \brief Generate and store all required local bonded interactions in \p idef and local exclusions in \p lexcls
 779  *
 780  * \returns Total count of bonded interactions in the local topology on this domain */
 781 static int make_local_bondeds_excls(gmx_domdec_t*           dd,
 782                                     gmx_domdec_zones_t*     zones,
 783                                     const gmx_mtop_t&       mtop,
 784                                     ArrayRef<const int64_t> atomInfo,
 785                                     const bool              checkDistanceMultiBody,
 786                                     const ivec              rcheck,
 787                                     const gmx_bool          checkDistanceTwoBody,
 788                                     const real              cutoff,
 789                                     const t_pbc*            pbc_null,
 790                                     ArrayRef<const RVec>    coordinates,
 791                                     InteractionDefinitions* idef,
 792                                     ListOfLists<int>*       lexcls)
 793 {
 794     int nzone_bondeds = 0;
 795
 796     if (dd->reverse_top->hasInterAtomicInteractions())
 797     {
 798         nzone_bondeds = zones->n;
 799     }
 800     else
 801     {
 802         /* Only single charge group (or atom) molecules, so interactions don't
 803          * cross zone boundaries and we only need to assign in the home zone.
 804          */
 805         nzone_bondeds = 1;
 806     }
 807
 808     /* We only use exclusions from i-zones to i- and j-zones */
 809     const int numIZonesForExclusions = (dd->haveExclusions ? zones->iZones.size() : 0);
 810
 811     gmx_reverse_top_t* rt = dd->reverse_top.get();
 812
 813     const real cutoffSquared = gmx::square(cutoff);
 814
 815     /* Clear the counts */
 816     idef->clear();
 817     int numBondedInteractions = 0;
 818
 819     lexcls->clear();
 820
 821     for (int izone = 0; izone < nzone_bondeds; izone++)
 822     {
 823         const int cg0 = zones->cg_range[izone];
 824         const int cg1 = zones->cg_range[izone + 1];
 825
 826         gmx::ArrayRef<thread_work_t> threadWorkObjects = rt->threadWorkObjects();
 827         const int                    numThreads        = threadWorkObjects.size();
 828 #pragma omp parallel for num_threads(numThreads) schedule(static)
 829         for (int thread = 0; thread < numThreads; thread++)
 830         {
 831             try
 832             {
 833                 InteractionDefinitions* idef_t = nullptr;
 834
 835                 int cg0t = cg0 + ((cg1 - cg0) * thread) / numThreads;
 836                 int cg1t = cg0 + ((cg1 - cg0) * (thread + 1)) / numThreads;
 837
 838                 if (thread == 0)
 839                 {
 840                     idef_t = idef;
 841                 }
 842                 else
 843                 {
 844                     idef_t = &threadWorkObjects[thread].idef;
 845                     idef_t->clear();
 846                 }
 847
 848                 threadWorkObjects[thread].numBondedInteractions =
 849                         make_bondeds_zone(rt,
 850                                           dd->globalAtomIndices,
 851                                           *dd->ga2la,
 852                                           zones,
 853                                           mtop.molblock,
 854                                           checkDistanceMultiBody,
 855                                           rcheck,
 856                                           checkDistanceTwoBody,
 857                                           cutoffSquared,
 858                                           pbc_null,
 859                                           coordinates,
 860                                           idef->iparams.data(),
 861                                           idef_t,
 862                                           izone,
 863                                           gmx::Range<int>(cg0t, cg1t));
 864
 865                 if (izone < numIZonesForExclusions)
 866                 {
 867                     ListOfLists<int>* excl_t = nullptr;
 868                     if (thread == 0)
 869                     {
 870                         // Thread 0 stores exclusions directly in the final storage
 871                         excl_t = lexcls;
 872                     }
 873                     else
 874                     {
 875                         // Threads > 0 store in temporary storage, starting at list index 0
 876                         excl_t = &threadWorkObjects[thread].excl;
 877                         excl_t->clear();
 878                     }
 879
 880                     /* No charge groups and no distance check required */
 881                     make_exclusions_zone(dd->globalAtomIndices,
 882                                          *dd->ga2la,
 883                                          zones,
 884                                          rt->molblockIndices(),
 885                                          mtop.moltype,
 886                                          atomInfo,
 887                                          excl_t,
 888                                          izone,
 889                                          cg0t,
 890                                          cg1t,
 891                                          mtop.intermolecularExclusionGroup);
 892                 }
 893             }
 894             GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
 895         }
 896
 897         if (threadWorkObjects.size() > 1)
 898         {
 899             combine_idef(idef, threadWorkObjects);
 900         }
 901
 902         for (const thread_work_t& th_work : threadWorkObjects)
 903         {
 904             numBondedInteractions += th_work.numBondedInteractions;
 905         }
 906
 907         if (izone < numIZonesForExclusions)
 908         {
 909             for (std::size_t th = 1; th < threadWorkObjects.size(); th++)
 910             {
 911                 lexcls->appendListOfLists(threadWorkObjects[th].excl);
 912             }
 913         }
 914     }
 915
 916     if (debug)
 917     {
 918         fprintf(debug, "We have %d exclusions\n", lexcls->numElements());
 919     }
 920
 921     return numBondedInteractions;
 922 }
 923
 924 int dd_make_local_top(gmx_domdec_t*        dd,
 925                       gmx_domdec_zones_t*  zones,
 926                       int                  npbcdim,
 927                       matrix               box,
 928                       rvec                 cellsize_min,
 929                       const ivec           npulse,
 930                       t_forcerec*          fr,
 931                       ArrayRef<const RVec> coordinates,
 932                       const gmx_mtop_t&    mtop,
 933                       gmx_localtop_t*      ltop)
 934 {
 935     real  rc = -1;
 936     ivec  rcheck;
 937     t_pbc pbc, *pbc_null = nullptr;
 938
 939     if (debug)
 940     {
 941         fprintf(debug, "Making local topology\n");
 942     }
 943
 944     bool checkDistanceMultiBody = false;
 945     bool checkDistanceTwoBody   = false;
 946
 947     if (dd->reverse_top->hasInterAtomicInteractions())
 948     {
 949         /* We need to check to which cell bondeds should be assigned */
 950         rc = dd_cutoff_twobody(dd);
 951         if (debug)
 952         {
 953             fprintf(debug, "Two-body bonded cut-off distance is %g\n", rc);
 954         }
 955
 956         /* Should we check distances when assigning bonded interactions? */
 957         for (int d = 0; d < DIM; d++)
 958         {
 959             rcheck[d] = FALSE;
 960             /* Only need to check for dimensions where the part of the box
 961              * that is not communicated is smaller than the cut-off.
 962              */
 963             if (d < npbcdim && dd->numCells[d] > 1
 964                 && (dd->numCells[d] - npulse[d]) * cellsize_min[d] < 2 * rc)
 965             {
 966                 if (dd->numCells[d] == 2)
 967                 {
 968                     rcheck[d]              = TRUE;
 969                     checkDistanceMultiBody = TRUE;
 970                 }
 971                 /* Check for interactions between two atoms,
 972                  * where we can allow interactions up to the cut-off,
 973                  * instead of up to the smallest cell dimension.
 974                  */
 975                 checkDistanceTwoBody = TRUE;
 976             }
 977             if (debug)
 978             {
 979                 fprintf(debug,
 980                         "dim %d cellmin %f bonded rcheck[%d] = %d, checkDistanceTwoBody = %s\n",
 981                         d,
 982                         cellsize_min[d],
 983                         d,
 984                         rcheck[d],
 985                         gmx::boolToString(checkDistanceTwoBody));
 986             }
 987         }
 988         if (checkDistanceMultiBody || checkDistanceTwoBody)
 989         {
 990             if (fr->bMolPBC)
 991             {
 992                 pbc_null = set_pbc_dd(&pbc, fr->pbcType, dd->numCells, TRUE, box);
 993             }
 994             else
 995             {
 996                 pbc_null = nullptr;
 997             }
 998         }
 999     }
1000
1001     int numBondedInteractionsToReduce = make_local_bondeds_excls(dd,
1002                                                                  zones,
1003                                                                  mtop,
1004                                                                  fr->atomInfo,
1005                                                                  checkDistanceMultiBody,
1006                                                                  rcheck,
1007                                                                  checkDistanceTwoBody,
1008                                                                  rc,
1009                                                                  pbc_null,
1010                                                                  coordinates,
1011                                                                  &ltop->idef,
1012                                                                  &ltop->excls);
1013
1014     /* The ilist is not sorted yet,
1015      * we can only do this when we have the charge arrays.
1016      */
1017     ltop->idef.ilsort = ilsortUNKNOWN;
1018
1019     return numBondedInteractionsToReduce;
1020 }
1021
1022 void dd_sort_local_top(const gmx_domdec_t& dd, gmx::ArrayRef<const int64_t> atomInfo, gmx_localtop_t* ltop)
1023 {
1024     if (dd.reverse_top->doSorting())
1025     {
1026         gmx_sort_ilist_fe(&ltop->idef, atomInfo);
1027     }
1028     else
1029     {
1030         ltop->idef.ilsort = ilsortNO_FE;
1031     }
1032 }