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