src/gromacs/domdec/domdec_topology.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, 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
  40  * while managing the construction, use and error checking for
  41  * topologies local to a DD rank.
  42  *
  43  * \author Berk Hess <hess@kth.se>
  44  * \ingroup module_domdec
  45  */
  46
  47 #include "gmxpre.h"
  48
  49 #include <cassert>
  50 #include <cstdlib>
  51 #include <cstring>
  52
  53 #include <algorithm>
  54 #include <memory>
  55 #include <string>
  56
  57 #include "gromacs/domdec/domdec.h"
  58 #include "gromacs/domdec/domdec_network.h"
  59 #include "gromacs/domdec/ga2la.h"
  60 #include "gromacs/gmxlib/network.h"
  61 #include "gromacs/math/vec.h"
  62 #include "gromacs/mdlib/forcerec.h"
  63 #include "gromacs/mdlib/gmx_omp_nthreads.h"
  64 #include "gromacs/mdlib/vsite.h"
  65 #include "gromacs/mdtypes/commrec.h"
  66 #include "gromacs/mdtypes/forcerec.h"
  67 #include "gromacs/mdtypes/inputrec.h"
  68 #include "gromacs/mdtypes/md_enums.h"
  69 #include "gromacs/mdtypes/mdatom.h"
  70 #include "gromacs/mdtypes/state.h"
  71 #include "gromacs/pbcutil/mshift.h"
  72 #include "gromacs/pbcutil/pbc.h"
  73 #include "gromacs/topology/ifunc.h"
  74 #include "gromacs/topology/mtop_util.h"
  75 #include "gromacs/topology/topsort.h"
  76 #include "gromacs/utility/cstringutil.h"
  77 #include "gromacs/utility/exceptions.h"
  78 #include "gromacs/utility/fatalerror.h"
  79 #include "gromacs/utility/gmxassert.h"
  80 #include "gromacs/utility/listoflists.h"
  81 #include "gromacs/utility/logger.h"
  82 #include "gromacs/utility/smalloc.h"
  83 #include "gromacs/utility/strconvert.h"
  84 #include "gromacs/utility/stringstream.h"
  85 #include "gromacs/utility/stringutil.h"
  86 #include "gromacs/utility/textwriter.h"
  87
  88 #include "domdec_constraints.h"
  89 #include "domdec_internal.h"
  90 #include "domdec_vsite.h"
  91 #include "dump.h"
  92
  93 using gmx::ArrayRef;
  94 using gmx::ListOfLists;
  95 using gmx::RVec;
  96
  97 /*! \brief The number of integer item in the local state, used for broadcasting of the state */
  98 #define NITEM_DD_INIT_LOCAL_STATE 5
  99
 100 struct reverse_ilist_t
 101 {
 102     std::vector<int> index;              /* Index for each atom into il          */
 103     std::vector<int> il;                 /* ftype|type|a0|...|an|ftype|...       */
 104     int              numAtomsInMolecule; /* The number of atoms in this molecule */
 105 };
 106
 107 struct MolblockIndices
 108 {
 109     int a_start;
 110     int a_end;
 111     int natoms_mol;
 112     int type;
 113 };
 114
 115 /*! \brief Struct for thread local work data for local topology generation */
 116 struct thread_work_t
 117 {
 118     /*! \brief Constructor
 119      *
 120      * \param[in] ffparams  The interaction parameters, the lifetime of the created object should not exceed the lifetime of the passed parameters
 121      */
 122     thread_work_t(const gmx_ffparams_t& ffparams) : idef(ffparams) {}
 123
 124     InteractionDefinitions         idef;               /**< Partial local topology */
 125     std::unique_ptr<gmx::VsitePbc> vsitePbc = nullptr; /**< vsite PBC structure */
 126     int                            nbonded  = 0;       /**< The number of bondeds in this struct */
 127     ListOfLists<int>               excl;               /**< List of exclusions */
 128     int                            excl_count = 0;     /**< The total exclusion count for \p excl */
 129 };
 130
 131 /*! \brief Struct for the reverse topology: links bonded interactions to atomsx */
 132 struct gmx_reverse_top_t
 133 {
 134     //! @cond Doxygen_Suppress
 135     //! \brief Are there constraints in this revserse top?
 136     bool bConstr = false;
 137     //! \brief Are there settles in this revserse top?
 138     bool bSettle = false;
 139     //! \brief All bonded interactions have to be assigned?
 140     bool bBCheck = false;
 141     //! \brief Are there bondeds/exclusions between atoms?
 142     bool bInterAtomicInteractions = false;
 143     //! \brief Reverse ilist for all moltypes
 144     std::vector<reverse_ilist_t> ril_mt;
 145     //! \brief The size of ril_mt[?].index summed over all entries
 146     int ril_mt_tot_size = 0;
 147     //! \brief The sorting state of bondeds for free energy
 148     int ilsort = ilsortUNKNOWN;
 149     //! \brief molblock to global atom index for quick lookup of molblocks on atom index
 150     std::vector<MolblockIndices> mbi;
 151
 152     //! \brief Do we have intermolecular interactions?
 153     bool bIntermolecularInteractions = false;
 154     //! \brief Intermolecular reverse ilist
 155     reverse_ilist_t ril_intermol;
 156
 157     /* Work data structures for multi-threading */
 158     //! \brief Thread work array for local topology generation
 159     std::vector<thread_work_t> th_work;
 160     //! @endcond
 161 };
 162
 163 /*! \brief Returns the number of atom entries for il in gmx_reverse_top_t */
 164 static int nral_rt(int ftype)
 165 {
 166     int nral;
 167
 168     nral = NRAL(ftype);
 169     if (interaction_function[ftype].flags & IF_VSITE)
 170     {
 171         /* With vsites the reverse topology contains an extra entry
 172          * for storing if constructing atoms are vsites.
 173          */
 174         nral += 1;
 175     }
 176
 177     return nral;
 178 }
 179
 180 /*! \brief Return whether interactions of type \p ftype need to be assigned exactly once */
 181 static gmx_bool dd_check_ftype(int ftype, gmx_bool bBCheck, gmx_bool bConstr, gmx_bool bSettle)
 182 {
 183     return ((((interaction_function[ftype].flags & IF_BOND) != 0U)
 184              && ((interaction_function[ftype].flags & IF_VSITE) == 0U)
 185              && (bBCheck || ((interaction_function[ftype].flags & IF_LIMZERO) == 0U)))
 186             || (bConstr && (ftype == F_CONSTR || ftype == F_CONSTRNC)) || (bSettle && ftype == F_SETTLE));
 187 }
 188
 189 /*! \brief Checks whether interactions have been assigned for one function type
 190  *
 191  * Loops over a list of interactions in the local topology of one function type
 192  * and flags each of the interactions as assigned in the global \p isAssigned list.
 193  * Exits with an inconsistency error when an interaction is assigned more than once.
 194  */
 195 static void flagInteractionsForType(const int                ftype,
 196                                     const InteractionList&   il,
 197                                     const reverse_ilist_t&   ril,
 198                                     const gmx::Range<int>&   atomRange,
 199                                     const int                numAtomsPerMolecule,
 200                                     gmx::ArrayRef<const int> globalAtomIndices,
 201                                     gmx::ArrayRef<int>       isAssigned)
 202 {
 203     const int nril_mol = ril.index[numAtomsPerMolecule];
 204     const int nral     = NRAL(ftype);
 205
 206     for (int i = 0; i < il.size(); i += 1 + nral)
 207     {
 208         // ia[0] is the interaction type, ia[1, ...] the atom indices
 209         const int* ia = il.iatoms.data() + i;
 210         // Extract the global atom index of the first atom in this interaction
 211         const int a0 = globalAtomIndices[ia[1]];
 212         /* Check if this interaction is in
 213          * the currently checked molblock.
 214          */
 215         if (atomRange.isInRange(a0))
 216         {
 217             // The molecule index in the list of this molecule type
 218             const int moleculeIndex = (a0 - atomRange.begin()) / numAtomsPerMolecule;
 219             const int atomOffset = (a0 - atomRange.begin()) - moleculeIndex * numAtomsPerMolecule;
 220             const int globalAtomStartInMolecule = atomRange.begin() + moleculeIndex * numAtomsPerMolecule;
 221             int       j_mol                     = ril.index[atomOffset];
 222             bool found                          = false;
 223             while (j_mol < ril.index[atomOffset + 1] && !found)
 224             {
 225                 const int j       = moleculeIndex * nril_mol + j_mol;
 226                 const int ftype_j = ril.il[j_mol];
 227                 /* Here we need to check if this interaction has
 228                  * not already been assigned, since we could have
 229                  * multiply defined interactions.
 230                  */
 231                 if (ftype == ftype_j && ia[0] == ril.il[j_mol + 1] && isAssigned[j] == 0)
 232                 {
 233                     /* Check the atoms */
 234                     found = true;
 235                     for (int a = 0; a < nral; a++)
 236                     {
 237                         if (globalAtomIndices[ia[1 + a]]
 238                             != globalAtomStartInMolecule + ril.il[j_mol + 2 + a])
 239                         {
 240                             found = false;
 241                         }
 242                     }
 243                     if (found)
 244                     {
 245                         isAssigned[j] = 1;
 246                     }
 247                 }
 248                 j_mol += 2 + nral_rt(ftype_j);
 249             }
 250             if (!found)
 251             {
 252                 gmx_incons("Some interactions seem to be assigned multiple times");
 253             }
 254         }
 255     }
 256 }
 257
 258 /*! \brief Help print error output when interactions are missing in a molblock
 259  *
 260  * \note This function needs to be called on all ranks (contains a global summation)
 261  */
 262 static std::string printMissingInteractionsMolblock(t_commrec*               cr,
 263                                                     const gmx_reverse_top_t& rt,
 264                                                     const char*              moltypename,
 265                                                     const reverse_ilist_t&   ril,
 266                                                     const gmx::Range<int>&   atomRange,
 267                                                     const int                numAtomsPerMolecule,
 268                                                     const int                numMolecules,
 269                                                     const InteractionDefinitions& idef)
 270 {
 271     const int               nril_mol = ril.index[numAtomsPerMolecule];
 272     std::vector<int>        isAssigned(numMolecules * nril_mol, 0);
 273     gmx::StringOutputStream stream;
 274     gmx::TextWriter         log(&stream);
 275
 276     for (int ftype = 0; ftype < F_NRE; ftype++)
 277     {
 278         if (dd_check_ftype(ftype, rt.bBCheck, rt.bConstr, rt.bSettle))
 279         {
 280             flagInteractionsForType(
 281                     ftype, idef.il[ftype], ril, atomRange, numAtomsPerMolecule, cr->dd->globalAtomIndices, isAssigned);
 282         }
 283     }
 284
 285     gmx_sumi(isAssigned.size(), isAssigned.data(), cr);
 286
 287     const int numMissingToPrint = 10;
 288     int       i                 = 0;
 289     for (int mol = 0; mol < numMolecules; mol++)
 290     {
 291         int j_mol = 0;
 292         while (j_mol < nril_mol)
 293         {
 294             int ftype = ril.il[j_mol];
 295             int nral  = NRAL(ftype);
 296             int j     = mol * nril_mol + j_mol;
 297             if (isAssigned[j] == 0 && !(interaction_function[ftype].flags & IF_VSITE))
 298             {
 299                 if (DDMASTER(cr->dd))
 300                 {
 301                     if (i == 0)
 302                     {
 303                         log.writeLineFormatted("Molecule type '%s'", moltypename);
 304                         log.writeLineFormatted(
 305                                 "the first %d missing interactions, except for exclusions:",
 306                                 numMissingToPrint);
 307                     }
 308                     log.writeStringFormatted("%20s atoms", interaction_function[ftype].longname);
 309                     int a;
 310                     for (a = 0; a < nral; a++)
 311                     {
 312                         log.writeStringFormatted("%5d", ril.il[j_mol + 2 + a] + 1);
 313                     }
 314                     while (a < 4)
 315                     {
 316                         log.writeString("     ");
 317                         a++;
 318                     }
 319                     log.writeString(" global");
 320                     for (a = 0; a < nral; a++)
 321                     {
 322                         log.writeStringFormatted("%6d",
 323                                                  atomRange.begin() + mol * numAtomsPerMolecule
 324                                                          + ril.il[j_mol + 2 + a] + 1);
 325                     }
 326                     log.ensureLineBreak();
 327                 }
 328                 i++;
 329                 if (i >= numMissingToPrint)
 330                 {
 331                     break;
 332                 }
 333             }
 334             j_mol += 2 + nral_rt(ftype);
 335         }
 336     }
 337
 338     return stream.toString();
 339 }
 340
 341 /*! \brief Help print error output when interactions are missing */
 342 static void printMissingInteractionsAtoms(const gmx::MDLogger&          mdlog,
 343                                           t_commrec*                    cr,
 344                                           const gmx_mtop_t&             mtop,
 345                                           const InteractionDefinitions& idef)
 346 {
 347     const gmx_reverse_top_t& rt = *cr->dd->reverse_top;
 348
 349     /* Print the atoms in the missing interactions per molblock */
 350     int a_end = 0;
 351     for (const gmx_molblock_t& molb : mtop.molblock)
 352     {
 353         const gmx_moltype_t& moltype = mtop.moltype[molb.type];
 354         const int            a_start = a_end;
 355         a_end                        = a_start + molb.nmol * moltype.atoms.nr;
 356         const gmx::Range<int> atomRange(a_start, a_end);
 357
 358         auto warning = printMissingInteractionsMolblock(
 359                 cr, rt, *(moltype.name), rt.ril_mt[molb.type], atomRange, moltype.atoms.nr, molb.nmol, idef);
 360
 361         GMX_LOG(mdlog.warning).appendText(warning);
 362     }
 363 }
 364
 365 void dd_print_missing_interactions(const gmx::MDLogger&           mdlog,
 366                                    t_commrec*                     cr,
 367                                    int                            local_count,
 368                                    const gmx_mtop_t*              top_global,
 369                                    const gmx_localtop_t*          top_local,
 370                                    gmx::ArrayRef<const gmx::RVec> x,
 371                                    const matrix                   box)
 372 {
 373     int           ndiff_tot, cl[F_NRE], n, ndiff, rest_global, rest_local;
 374     int           ftype, nral;
 375     gmx_domdec_t* dd;
 376
 377     dd = cr->dd;
 378
 379     GMX_LOG(mdlog.warning)
 380             .appendText(
 381                     "Not all bonded interactions have been properly assigned to the domain "
 382                     "decomposition cells");
 383
 384     ndiff_tot = local_count - dd->nbonded_global;
 385
 386     for (ftype = 0; ftype < F_NRE; ftype++)
 387     {
 388         nral      = NRAL(ftype);
 389         cl[ftype] = top_local->idef.il[ftype].size() / (1 + nral);
 390     }
 391
 392     gmx_sumi(F_NRE, cl, cr);
 393
 394     if (DDMASTER(dd))
 395     {
 396         GMX_LOG(mdlog.warning).appendText("A list of missing interactions:");
 397         rest_global = dd->nbonded_global;
 398         rest_local  = local_count;
 399         for (ftype = 0; ftype < F_NRE; ftype++)
 400         {
 401             /* In the reverse and local top all constraints are merged
 402              * into F_CONSTR. So in the if statement we skip F_CONSTRNC
 403              * and add these constraints when doing F_CONSTR.
 404              */
 405             if (((interaction_function[ftype].flags & IF_BOND)
 406                  && (dd->reverse_top->bBCheck || !(interaction_function[ftype].flags & IF_LIMZERO)))
 407                 || (dd->reverse_top->bConstr && ftype == F_CONSTR)
 408                 || (dd->reverse_top->bSettle && ftype == F_SETTLE))
 409             {
 410                 n = gmx_mtop_ftype_count(top_global, ftype);
 411                 if (ftype == F_CONSTR)
 412                 {
 413                     n += gmx_mtop_ftype_count(top_global, F_CONSTRNC);
 414                 }
 415                 ndiff = cl[ftype] - n;
 416                 if (ndiff != 0)
 417                 {
 418                     GMX_LOG(mdlog.warning)
 419                             .appendTextFormatted("%20s of %6d missing %6d",
 420                                                  interaction_function[ftype].longname,
 421                                                  n,
 422                                                  -ndiff);
 423                 }
 424                 rest_global -= n;
 425                 rest_local -= cl[ftype];
 426             }
 427         }
 428
 429         ndiff = rest_local - rest_global;
 430         if (ndiff != 0)
 431         {
 432             GMX_LOG(mdlog.warning).appendTextFormatted("%20s of %6d missing %6d", "exclusions", rest_global, -ndiff);
 433         }
 434     }
 435
 436     printMissingInteractionsAtoms(mdlog, cr, *top_global, top_local->idef);
 437     write_dd_pdb("dd_dump_err", 0, "dump", top_global, cr, -1, as_rvec_array(x.data()), box);
 438
 439     std::string errorMessage;
 440
 441     if (ndiff_tot > 0)
 442     {
 443         errorMessage =
 444                 "One or more interactions were assigned to multiple domains of the domain "
 445                 "decompostion. Please report this bug.";
 446     }
 447     else
 448     {
 449         errorMessage = gmx::formatString(
 450                 "%d of the %d bonded interactions could not be calculated because some atoms "
 451                 "involved moved further apart than the multi-body cut-off distance (%g nm) or the "
 452                 "two-body cut-off distance (%g nm), see option -rdd, for pairs and tabulated bonds "
 453                 "also see option -ddcheck",
 454                 -ndiff_tot,
 455                 cr->dd->nbonded_global,
 456                 dd_cutoff_multibody(dd),
 457                 dd_cutoff_twobody(dd));
 458     }
 459     gmx_fatal_collective(FARGS, cr->mpi_comm_mygroup, MASTER(cr), "%s", errorMessage.c_str());
 460 }
 461
 462 /*! \brief Return global topology molecule information for global atom index \p i_gl */
 463 static void global_atomnr_to_moltype_ind(const gmx_reverse_top_t* rt, int i_gl, int* mb, int* mt, int* mol, int* i_mol)
 464 {
 465     const MolblockIndices* mbi   = rt->mbi.data();
 466     int                    start = 0;
 467     int                    end   = rt->mbi.size(); /* exclusive */
 468     int                    mid;
 469
 470     /* binary search for molblock_ind */
 471     while (TRUE)
 472     {
 473         mid = (start + end) / 2;
 474         if (i_gl >= mbi[mid].a_end)
 475         {
 476             start = mid + 1;
 477         }
 478         else if (i_gl < mbi[mid].a_start)
 479         {
 480             end = mid;
 481         }
 482         else
 483         {
 484             break;
 485         }
 486     }
 487
 488     *mb = mid;
 489     mbi += mid;
 490
 491     *mt    = mbi->type;
 492     *mol   = (i_gl - mbi->a_start) / mbi->natoms_mol;
 493     *i_mol = (i_gl - mbi->a_start) - (*mol) * mbi->natoms_mol;
 494 }
 495
 496 /*! \brief Returns the maximum number of exclusions per atom */
 497 static int getMaxNumExclusionsPerAtom(const ListOfLists<int>& excls)
 498 {
 499     int maxNumExcls = 0;
 500     for (gmx::index at = 0; at < excls.ssize(); at++)
 501     {
 502         const auto list     = excls[at];
 503         const int  numExcls = list.ssize();
 504
 505         GMX_RELEASE_ASSERT(numExcls != 1 || list[0] == at,
 506                            "With 1 exclusion we expect a self-exclusion");
 507
 508         maxNumExcls = std::max(maxNumExcls, numExcls);
 509     }
 510
 511     return maxNumExcls;
 512 }
 513
 514 /*! \brief Run the reverse ilist generation and store it in r_il when \p bAssign = TRUE */
 515 static int low_make_reverse_ilist(const InteractionLists&  il_mt,
 516                                   const t_atom*            atom,
 517                                   int*                     count,
 518                                   gmx_bool                 bConstr,
 519                                   gmx_bool                 bSettle,
 520                                   gmx_bool                 bBCheck,
 521                                   gmx::ArrayRef<const int> r_index,
 522                                   gmx::ArrayRef<int>       r_il,
 523                                   gmx_bool                 bLinkToAllAtoms,
 524                                   gmx_bool                 bAssign)
 525 {
 526     int ftype, j, nlink, link;
 527     int a;
 528     int nint;
 529
 530     nint = 0;
 531     for (ftype = 0; ftype < F_NRE; ftype++)
 532     {
 533         if ((interaction_function[ftype].flags & (IF_BOND | IF_VSITE))
 534             || (bConstr && (ftype == F_CONSTR || ftype == F_CONSTRNC)) || (bSettle && ftype == F_SETTLE))
 535         {
 536             const bool  bVSite = ((interaction_function[ftype].flags & IF_VSITE) != 0U);
 537             const int   nral   = NRAL(ftype);
 538             const auto& il     = il_mt[ftype];
 539             for (int i = 0; i < il.size(); i += 1 + nral)
 540             {
 541                 const int* ia = il.iatoms.data() + i;
 542                 if (bLinkToAllAtoms)
 543                 {
 544                     if (bVSite)
 545                     {
 546                         /* We don't need the virtual sites for the cg-links */
 547                         nlink = 0;
 548                     }
 549                     else
 550                     {
 551                         nlink = nral;
 552                     }
 553                 }
 554                 else
 555                 {
 556                     /* Couple to the first atom in the interaction */
 557                     nlink = 1;
 558                 }
 559                 for (link = 0; link < nlink; link++)
 560                 {
 561                     a = ia[1 + link];
 562                     if (bAssign)
 563                     {
 564                         GMX_ASSERT(!r_il.empty(), "with bAssign not allowed to be empty");
 565                         GMX_ASSERT(!r_index.empty(), "with bAssign not allowed to be empty");
 566                         r_il[r_index[a] + count[a]]     = (ftype == F_CONSTRNC ? F_CONSTR : ftype);
 567                         r_il[r_index[a] + count[a] + 1] = ia[0];
 568                         for (j = 1; j < 1 + nral; j++)
 569                         {
 570                             /* Store the molecular atom number */
 571                             r_il[r_index[a] + count[a] + 1 + j] = ia[j];
 572                         }
 573                     }
 574                     if (interaction_function[ftype].flags & IF_VSITE)
 575                     {
 576                         if (bAssign)
 577                         {
 578                             /* Add an entry to iatoms for storing
 579                              * which of the constructing atoms are
 580                              * vsites again.
 581                              */
 582                             r_il[r_index[a] + count[a] + 2 + nral] = 0;
 583                             for (j = 2; j < 1 + nral; j++)
 584                             {
 585                                 if (atom[ia[j]].ptype == eptVSite)
 586                                 {
 587                                     r_il[r_index[a] + count[a] + 2 + nral] |= (2 << j);
 588                                 }
 589                             }
 590                         }
 591                     }
 592                     else
 593                     {
 594                         /* We do not count vsites since they are always
 595                          * uniquely assigned and can be assigned
 596                          * to multiple nodes with recursive vsites.
 597                          */
 598                         if (bBCheck || !(interaction_function[ftype].flags & IF_LIMZERO))
 599                         {
 600                             nint++;
 601                         }
 602                     }
 603                     count[a] += 2 + nral_rt(ftype);
 604                 }
 605             }
 606         }
 607     }
 608
 609     return nint;
 610 }
 611
 612 /*! \brief Make the reverse ilist: a list of bonded interactions linked to atoms */
 613 static int make_reverse_ilist(const InteractionLists& ilist,
 614                               const t_atoms*          atoms,
 615                               gmx_bool                bConstr,
 616                               gmx_bool                bSettle,
 617                               gmx_bool                bBCheck,
 618                               gmx_bool                bLinkToAllAtoms,
 619                               reverse_ilist_t*        ril_mt)
 620 {
 621     int nat_mt, *count, i, nint_mt;
 622
 623     /* Count the interactions */
 624     nat_mt = atoms->nr;
 625     snew(count, nat_mt);
 626     low_make_reverse_ilist(
 627             ilist, atoms->atom, count, bConstr, bSettle, bBCheck, {}, {}, bLinkToAllAtoms, FALSE);
 628
 629     ril_mt->index.push_back(0);
 630     for (i = 0; i < nat_mt; i++)
 631     {
 632         ril_mt->index.push_back(ril_mt->index[i] + count[i]);
 633         count[i] = 0;
 634     }
 635     ril_mt->il.resize(ril_mt->index[nat_mt]);
 636
 637     /* Store the interactions */
 638     nint_mt = low_make_reverse_ilist(
 639             ilist, atoms->atom, count, bConstr, bSettle, bBCheck, ril_mt->index, ril_mt->il, bLinkToAllAtoms, TRUE);
 640
 641     sfree(count);
 642
 643     ril_mt->numAtomsInMolecule = atoms->nr;
 644
 645     return nint_mt;
 646 }
 647
 648 /*! \brief Generate the reverse topology */
 649 static gmx_reverse_top_t make_reverse_top(const gmx_mtop_t* mtop,
 650                                           gmx_bool          bFE,
 651                                           gmx_bool          bConstr,
 652                                           gmx_bool          bSettle,
 653                                           gmx_bool          bBCheck,
 654                                           int*              nint)
 655 {
 656     gmx_reverse_top_t rt;
 657
 658     /* Should we include constraints (for SHAKE) in rt? */
 659     rt.bConstr = bConstr;
 660     rt.bSettle = bSettle;
 661     rt.bBCheck = bBCheck;
 662
 663     rt.bInterAtomicInteractions = mtop->bIntermolecularInteractions;
 664     rt.ril_mt.resize(mtop->moltype.size());
 665     rt.ril_mt_tot_size = 0;
 666     std::vector<int> nint_mt;
 667     for (size_t mt = 0; mt < mtop->moltype.size(); mt++)
 668     {
 669         const gmx_moltype_t& molt = mtop->moltype[mt];
 670         if (molt.atoms.nr > 1)
 671         {
 672             rt.bInterAtomicInteractions = true;
 673         }
 674
 675         /* Make the atom to interaction list for this molecule type */
 676         int numberOfInteractions = make_reverse_ilist(
 677                 molt.ilist, &molt.atoms, rt.bConstr, rt.bSettle, rt.bBCheck, FALSE, &rt.ril_mt[mt]);
 678         nint_mt.push_back(numberOfInteractions);
 679
 680         rt.ril_mt_tot_size += rt.ril_mt[mt].index[molt.atoms.nr];
 681     }
 682     if (debug)
 683     {
 684         fprintf(debug, "The total size of the atom to interaction index is %d integers\n", rt.ril_mt_tot_size);
 685     }
 686
 687     *nint = 0;
 688     for (const gmx_molblock_t& molblock : mtop->molblock)
 689     {
 690         *nint += molblock.nmol * nint_mt[molblock.type];
 691     }
 692
 693     /* Make an intermolecular reverse top, if necessary */
 694     rt.bIntermolecularInteractions = mtop->bIntermolecularInteractions;
 695     if (rt.bIntermolecularInteractions)
 696     {
 697         t_atoms atoms_global;
 698
 699         atoms_global.nr   = mtop->natoms;
 700         atoms_global.atom = nullptr; /* Only used with virtual sites */
 701
 702         GMX_RELEASE_ASSERT(mtop->intermolecular_ilist,
 703                            "We should have an ilist when intermolecular interactions are on");
 704
 705         *nint += make_reverse_ilist(
 706                 *mtop->intermolecular_ilist, &atoms_global, rt.bConstr, rt.bSettle, rt.bBCheck, FALSE, &rt.ril_intermol);
 707     }
 708
 709     if (bFE && gmx_mtop_bondeds_free_energy(mtop))
 710     {
 711         rt.ilsort = ilsortFE_UNSORTED;
 712     }
 713     else
 714     {
 715         rt.ilsort = ilsortNO_FE;
 716     }
 717
 718     /* Make a molblock index for fast searching */
 719     int i = 0;
 720     for (size_t mb = 0; mb < mtop->molblock.size(); mb++)
 721     {
 722         const gmx_molblock_t& molb           = mtop->molblock[mb];
 723         const int             numAtomsPerMol = mtop->moltype[molb.type].atoms.nr;
 724         MolblockIndices       mbi;
 725         mbi.a_start = i;
 726         i += molb.nmol * numAtomsPerMol;
 727         mbi.a_end      = i;
 728         mbi.natoms_mol = numAtomsPerMol;
 729         mbi.type       = molb.type;
 730         rt.mbi.push_back(mbi);
 731     }
 732
 733     for (int th = 0; th < gmx_omp_nthreads_get(emntDomdec); th++)
 734     {
 735         rt.th_work.emplace_back(mtop->ffparams);
 736     }
 737
 738     return rt;
 739 }
 740
 741 void dd_make_reverse_top(FILE*                           fplog,
 742                          gmx_domdec_t*                   dd,
 743                          const gmx_mtop_t*               mtop,
 744                          const gmx::VirtualSitesHandler* vsite,
 745                          const t_inputrec*               ir,
 746                          gmx_bool                        bBCheck)
 747 {
 748     if (fplog)
 749     {
 750         fprintf(fplog, "\nLinking all bonded interactions to atoms\n");
 751     }
 752
 753     /* If normal and/or settle constraints act only within charge groups,
 754      * we can store them in the reverse top and simply assign them to domains.
 755      * Otherwise we need to assign them to multiple domains and set up
 756      * the parallel version constraint algorithm(s).
 757      */
 758
 759     dd->reverse_top  = new gmx_reverse_top_t;
 760     *dd->reverse_top = make_reverse_top(mtop,
 761                                         ir->efep != efepNO,
 762                                         !dd->comm->systemInfo.haveSplitConstraints,
 763                                         !dd->comm->systemInfo.haveSplitSettles,
 764                                         bBCheck,
 765                                         &dd->nbonded_global);
 766
 767     dd->haveExclusions = false;
 768     for (const gmx_molblock_t& molb : mtop->molblock)
 769     {
 770         const int maxNumExclusionsPerAtom = getMaxNumExclusionsPerAtom(mtop->moltype[molb.type].excls);
 771         // We checked above that max 1 exclusion means only self exclusions
 772         if (maxNumExclusionsPerAtom > 1)
 773         {
 774             dd->haveExclusions = true;
 775         }
 776     }
 777
 778     if (vsite && vsite->numInterUpdategroupVirtualSites() > 0)
 779     {
 780         if (fplog)
 781         {
 782             fprintf(fplog,
 783                     "There are %d inter update-group virtual sites,\n"
 784                     "will an extra communication step for selected coordinates and forces\n",
 785                     vsite->numInterUpdategroupVirtualSites());
 786         }
 787         init_domdec_vsites(dd, vsite->numInterUpdategroupVirtualSites());
 788     }
 789
 790     if (dd->comm->systemInfo.haveSplitConstraints || dd->comm->systemInfo.haveSplitSettles)
 791     {
 792         init_domdec_constraints(dd, mtop);
 793     }
 794     if (fplog)
 795     {
 796         fprintf(fplog, "\n");
 797     }
 798 }
 799
 800 /*! \brief Store a vsite interaction at the end of \p il
 801  *
 802  * This routine is very similar to add_ifunc, but vsites interactions
 803  * have more work to do than other kinds of interactions, and the
 804  * complex way nral (and thus vector contents) depends on ftype
 805  * confuses static analysis tools unless we fuse the vsite
 806  * atom-indexing organization code with the ifunc-adding code, so that
 807  * they can see that nral is the same value. */
 808 static inline void add_ifunc_for_vsites(t_iatom*           tiatoms,
 809                                         const gmx_ga2la_t& ga2la,
 810                                         int                nral,
 811                                         gmx_bool           bHomeA,
 812                                         int                a,
 813                                         int                a_gl,
 814                                         int                a_mol,
 815                                         const t_iatom*     iatoms,
 816                                         InteractionList*   il)
 817 {
 818     /* Copy the type */
 819     tiatoms[0] = iatoms[0];
 820
 821     if (bHomeA)
 822     {
 823         /* We know the local index of the first atom */
 824         tiatoms[1] = a;
 825     }
 826     else
 827     {
 828         /* Convert later in make_local_vsites */
 829         tiatoms[1] = -a_gl - 1;
 830     }
 831
 832     GMX_ASSERT(nral >= 2 && nral <= 5, "Invalid nral for vsites");
 833     for (int k = 2; k < 1 + nral; k++)
 834     {
 835         int ak_gl = a_gl + iatoms[k] - a_mol;
 836         if (const int* homeIndex = ga2la.findHome(ak_gl))
 837         {
 838             tiatoms[k] = *homeIndex;
 839         }
 840         else
 841         {
 842             /* Copy the global index, convert later in make_local_vsites */
 843             tiatoms[k] = -(ak_gl + 1);
 844         }
 845         // Note that ga2la_get_home always sets the third parameter if
 846         // it returns TRUE
 847     }
 848     il->push_back(tiatoms[0], nral, tiatoms + 1);
 849 }
 850
 851 /*! \brief Store a position restraint in idef and iatoms, complex because the parameters are different for each entry */
 852 static void add_posres(int                     mol,
 853                        int                     a_mol,
 854                        int                     numAtomsInMolecule,
 855                        const gmx_molblock_t*   molb,
 856                        t_iatom*                iatoms,
 857                        const t_iparams*        ip_in,
 858                        InteractionDefinitions* idef)
 859 {
 860     /* This position restraint has not been added yet,
 861      * so it's index is the current number of position restraints.
 862      */
 863     const int n = idef->il[F_POSRES].size() / 2;
 864
 865     /* Get the position restraint coordinates from the molblock */
 866     const int a_molb = mol * numAtomsInMolecule + a_mol;
 867     GMX_ASSERT(a_molb < ssize(molb->posres_xA),
 868                "We need a sufficient number of position restraint coordinates");
 869     /* Copy the force constants */
 870     t_iparams ip        = ip_in[iatoms[0]];
 871     ip.posres.pos0A[XX] = molb->posres_xA[a_molb][XX];
 872     ip.posres.pos0A[YY] = molb->posres_xA[a_molb][YY];
 873     ip.posres.pos0A[ZZ] = molb->posres_xA[a_molb][ZZ];
 874     if (!molb->posres_xB.empty())
 875     {
 876         ip.posres.pos0B[XX] = molb->posres_xB[a_molb][XX];
 877         ip.posres.pos0B[YY] = molb->posres_xB[a_molb][YY];
 878         ip.posres.pos0B[ZZ] = molb->posres_xB[a_molb][ZZ];
 879     }
 880     else
 881     {
 882         ip.posres.pos0B[XX] = ip.posres.pos0A[XX];
 883         ip.posres.pos0B[YY] = ip.posres.pos0A[YY];
 884         ip.posres.pos0B[ZZ] = ip.posres.pos0A[ZZ];
 885     }
 886     /* Set the parameter index for idef->iparams_posres */
 887     iatoms[0] = n;
 888     idef->iparams_posres.push_back(ip);
 889
 890     GMX_ASSERT(int(idef->iparams_posres.size()) == n + 1,
 891                "The index of the parameter type should match n");
 892 }
 893
 894 /*! \brief Store a flat-bottomed position restraint in idef and iatoms, complex because the parameters are different for each entry */
 895 static void add_fbposres(int                     mol,
 896                          int                     a_mol,
 897                          int                     numAtomsInMolecule,
 898                          const gmx_molblock_t*   molb,
 899                          t_iatom*                iatoms,
 900                          const t_iparams*        ip_in,
 901                          InteractionDefinitions* idef)
 902 {
 903     /* This flat-bottom position restraint has not been added yet,
 904      * so it's index is the current number of position restraints.
 905      */
 906     const int n = idef->il[F_FBPOSRES].size() / 2;
 907
 908     /* Get the position restraint coordinats from the molblock */
 909     const int a_molb = mol * numAtomsInMolecule + a_mol;
 910     GMX_ASSERT(a_molb < ssize(molb->posres_xA),
 911                "We need a sufficient number of position restraint coordinates");
 912     /* Copy the force constants */
 913     t_iparams ip = ip_in[iatoms[0]];
 914     /* Take reference positions from A position of normal posres */
 915     ip.fbposres.pos0[XX] = molb->posres_xA[a_molb][XX];
 916     ip.fbposres.pos0[YY] = molb->posres_xA[a_molb][YY];
 917     ip.fbposres.pos0[ZZ] = molb->posres_xA[a_molb][ZZ];
 918
 919     /* Note: no B-type for flat-bottom posres */
 920
 921     /* Set the parameter index for idef->iparams_fbposres */
 922     iatoms[0] = n;
 923     idef->iparams_fbposres.push_back(ip);
 924
 925     GMX_ASSERT(int(idef->iparams_fbposres.size()) == n + 1,
 926                "The index of the parameter type should match n");
 927 }
 928
 929 /*! \brief Store a virtual site interaction, complex because of PBC and recursion */
 930 static void add_vsite(const gmx_ga2la_t&       ga2la,
 931                       gmx::ArrayRef<const int> index,
 932                       gmx::ArrayRef<const int> rtil,
 933                       int                      ftype,
 934                       int                      nral,
 935                       gmx_bool                 bHomeA,
 936                       int                      a,
 937                       int                      a_gl,
 938                       int                      a_mol,
 939                       const t_iatom*           iatoms,
 940                       InteractionDefinitions*  idef)
 941 {
 942     int     k;
 943     t_iatom tiatoms[1 + MAXATOMLIST];
 944     int     j, ftype_r, nral_r;
 945
 946     /* Add this interaction to the local topology */
 947     add_ifunc_for_vsites(tiatoms, ga2la, nral, bHomeA, a, a_gl, a_mol, iatoms, &idef->il[ftype]);
 948
 949     if (iatoms[1 + nral])
 950     {
 951         /* Check for recursion */
 952         for (k = 2; k < 1 + nral; k++)
 953         {
 954             if ((iatoms[1 + nral] & (2 << k)) && (tiatoms[k] < 0))
 955             {
 956                 /* This construction atoms is a vsite and not a home atom */
 957                 if (gmx_debug_at)
 958                 {
 959                     fprintf(debug,
 960                             "Constructing atom %d of vsite atom %d is a vsite and non-home\n",
 961                             iatoms[k] + 1,
 962                             a_mol + 1);
 963                 }
 964                 /* Find the vsite construction */
 965
 966                 /* Check all interactions assigned to this atom */
 967                 j = index[iatoms[k]];
 968                 while (j < index[iatoms[k] + 1])
 969                 {
 970                     ftype_r = rtil[j++];
 971                     nral_r  = NRAL(ftype_r);
 972                     if (interaction_function[ftype_r].flags & IF_VSITE)
 973                     {
 974                         /* Add this vsite (recursion) */
 975                         add_vsite(ga2la,
 976                                   index,
 977                                   rtil,
 978                                   ftype_r,
 979                                   nral_r,
 980                                   FALSE,
 981                                   -1,
 982                                   a_gl + iatoms[k] - iatoms[1],
 983                                   iatoms[k],
 984                                   rtil.data() + j,
 985                                   idef);
 986                     }
 987                     j += 1 + nral_rt(ftype_r);
 988                 }
 989             }
 990         }
 991     }
 992 }
 993
 994 /*! \brief Returns the squared distance between atoms \p i and \p j */
 995 static real dd_dist2(t_pbc* pbc_null, const rvec* x, const int i, int j)
 996 {
 997     rvec dx;
 998
 999     if (pbc_null)
1000     {
1001         pbc_dx_aiuc(pbc_null, x[i], x[j], dx);
1002     }
1003     else
1004     {
1005         rvec_sub(x[i], x[j], dx);
1006     }
1007
1008     return norm2(dx);
1009 }
1010
1011 /*! \brief Append t_idef structures 1 to nsrc in src to *dest */
1012 static void combine_idef(InteractionDefinitions* dest, gmx::ArrayRef<const thread_work_t> src)
1013 {
1014     int ftype;
1015
1016     for (ftype = 0; ftype < F_NRE; ftype++)
1017     {
1018         int n = 0;
1019         for (gmx::index s = 1; s < src.ssize(); s++)
1020         {
1021             n += src[s].idef.il[ftype].size();
1022         }
1023         if (n > 0)
1024         {
1025             for (gmx::index s = 1; s < src.ssize(); s++)
1026             {
1027                 dest->il[ftype].append(src[s].idef.il[ftype]);
1028             }
1029
1030             /* Position restraints need an additional treatment */
1031             if (ftype == F_POSRES || ftype == F_FBPOSRES)
1032             {
1033                 int                     nposres = dest->il[ftype].size() / 2;
1034                 std::vector<t_iparams>& iparams_dest =
1035                         (ftype == F_POSRES ? dest->iparams_posres : dest->iparams_fbposres);
1036
1037                 /* Set nposres to the number of original position restraints in dest */
1038                 for (gmx::index s = 1; s < src.ssize(); s++)
1039                 {
1040                     nposres -= src[s].idef.il[ftype].size() / 2;
1041                 }
1042
1043                 for (gmx::index s = 1; s < src.ssize(); s++)
1044                 {
1045                     const std::vector<t_iparams>& iparams_src =
1046                             (ftype == F_POSRES ? src[s].idef.iparams_posres : src[s].idef.iparams_fbposres);
1047                     iparams_dest.insert(iparams_dest.end(), iparams_src.begin(), iparams_src.end());
1048
1049                     /* Correct the indices into iparams_posres */
1050                     for (int i = 0; i < src[s].idef.il[ftype].size() / 2; i++)
1051                     {
1052                         /* Correct the index into iparams_posres */
1053                         dest->il[ftype].iatoms[nposres * 2] = nposres;
1054                         nposres++;
1055                     }
1056                 }
1057                 GMX_RELEASE_ASSERT(
1058                         int(iparams_dest.size()) == nposres,
1059                         "The number of parameters should match the number of restraints");
1060             }
1061         }
1062     }
1063 }
1064
1065 /*! \brief Check and when available assign bonded interactions for local atom i
1066  */
1067 static inline void check_assign_interactions_atom(int                       i,
1068                                                   int                       i_gl,
1069                                                   int                       mol,
1070                                                   int                       i_mol,
1071                                                   int                       numAtomsInMolecule,
1072                                                   gmx::ArrayRef<const int>  index,
1073                                                   gmx::ArrayRef<const int>  rtil,
1074                                                   gmx_bool                  bInterMolInteractions,
1075                                                   int                       ind_start,
1076                                                   int                       ind_end,
1077                                                   const gmx_domdec_t*       dd,
1078                                                   const gmx_domdec_zones_t* zones,
1079                                                   const gmx_molblock_t*     molb,
1080                                                   gmx_bool                  bRCheckMB,
1081                                                   const ivec                rcheck,
1082                                                   gmx_bool                  bRCheck2B,
1083                                                   real                      rc2,
1084                                                   t_pbc*                    pbc_null,
1085                                                   rvec*                     cg_cm,
1086                                                   const t_iparams*          ip_in,
1087                                                   InteractionDefinitions*   idef,
1088                                                   int                       iz,
1089                                                   gmx_bool                  bBCheck,
1090                                                   int*                      nbonded_local)
1091 {
1092     gmx::ArrayRef<const DDPairInteractionRanges> iZones = zones->iZones;
1093
1094     int j = ind_start;
1095     while (j < ind_end)
1096     {
1097         t_iatom tiatoms[1 + MAXATOMLIST];
1098
1099         const int ftype  = rtil[j++];
1100         auto      iatoms = gmx::constArrayRefFromArray(rtil.data() + j, rtil.size() - j);
1101         const int nral   = NRAL(ftype);
1102         if (interaction_function[ftype].flags & IF_VSITE)
1103         {
1104             assert(!bInterMolInteractions);
1105             /* The vsite construction goes where the vsite itself is */
1106             if (iz == 0)
1107             {
1108                 add_vsite(*dd->ga2la, index, rtil, ftype, nral, TRUE, i, i_gl, i_mol, iatoms.data(), idef);
1109             }
1110         }
1111         else
1112         {
1113             gmx_bool bUse;
1114
1115             /* Copy the type */
1116             tiatoms[0] = iatoms[0];
1117
1118             if (nral == 1)
1119             {
1120                 assert(!bInterMolInteractions);
1121                 /* Assign single-body interactions to the home zone */
1122                 if (iz == 0)
1123                 {
1124                     bUse       = TRUE;
1125                     tiatoms[1] = i;
1126                     if (ftype == F_POSRES)
1127                     {
1128                         add_posres(mol, i_mol, numAtomsInMolecule, molb, tiatoms, ip_in, idef);
1129                     }
1130                     else if (ftype == F_FBPOSRES)
1131                     {
1132                         add_fbposres(mol, i_mol, numAtomsInMolecule, molb, tiatoms, ip_in, idef);
1133                     }
1134                 }
1135                 else
1136                 {
1137                     bUse = FALSE;
1138                 }
1139             }
1140             else if (nral == 2)
1141             {
1142                 /* This is a two-body interaction, we can assign
1143                  * analogous to the non-bonded assignments.
1144                  */
1145                 int k_gl;
1146
1147                 if (!bInterMolInteractions)
1148                 {
1149                     /* Get the global index using the offset in the molecule */
1150                     k_gl = i_gl + iatoms[2] - i_mol;
1151                 }
1152                 else
1153                 {
1154                     k_gl = iatoms[2];
1155                 }
1156                 if (const auto* entry = dd->ga2la->find(k_gl))
1157                 {
1158                     int kz = entry->cell;
1159                     if (kz >= zones->n)
1160                     {
1161                         kz -= zones->n;
1162                     }
1163                     /* Check zone interaction assignments */
1164                     bUse = ((iz < iZones.ssize() && iz <= kz && iZones[iz].jZoneRange.isInRange(kz))
1165                             || (kz < iZones.ssize() && iz > kz && iZones[kz].jZoneRange.isInRange(iz)));
1166                     if (bUse)
1167                     {
1168                         GMX_ASSERT(ftype != F_CONSTR || (iz == 0 && kz == 0),
1169                                    "Constraint assigned here should only involve home atoms");
1170
1171                         tiatoms[1] = i;
1172                         tiatoms[2] = entry->la;
1173                         /* If necessary check the cgcm distance */
1174                         if (bRCheck2B && dd_dist2(pbc_null, cg_cm, tiatoms[1], tiatoms[2]) >= rc2)
1175                         {
1176                             bUse = FALSE;
1177                         }
1178                     }
1179                 }
1180                 else
1181                 {
1182                     bUse = false;
1183                 }
1184             }
1185             else
1186             {
1187                 /* Assign this multi-body bonded interaction to
1188                  * the local node if we have all the atoms involved
1189                  * (local or communicated) and the minimum zone shift
1190                  * in each dimension is zero, for dimensions
1191                  * with 2 DD cells an extra check may be necessary.
1192                  */
1193                 ivec k_zero, k_plus;
1194                 int  k;
1195
1196                 bUse = TRUE;
1197                 clear_ivec(k_zero);
1198                 clear_ivec(k_plus);
1199                 for (k = 1; k <= nral && bUse; k++)
1200                 {
1201                     int k_gl;
1202                     if (!bInterMolInteractions)
1203                     {
1204                         /* Get the global index using the offset in the molecule */
1205                         k_gl = i_gl + iatoms[k] - i_mol;
1206                     }
1207                     else
1208                     {
1209                         k_gl = iatoms[k];
1210                     }
1211                     const auto* entry = dd->ga2la->find(k_gl);
1212                     if (entry == nullptr || entry->cell >= zones->n)
1213                     {
1214                         /* We do not have this atom of this interaction
1215                          * locally, or it comes from more than one cell
1216                          * away.
1217                          */
1218                         bUse = FALSE;
1219                     }
1220                     else
1221                     {
1222                         int d;
1223
1224                         tiatoms[k] = entry->la;
1225                         for (d = 0; d < DIM; d++)
1226                         {
1227                             if (zones->shift[entry->cell][d] == 0)
1228                             {
1229                                 k_zero[d] = k;
1230                             }
1231                             else
1232                             {
1233                                 k_plus[d] = k;
1234                             }
1235                         }
1236                     }
1237                 }
1238                 bUse = (bUse && (k_zero[XX] != 0) && (k_zero[YY] != 0) && (k_zero[ZZ] != 0));
1239                 if (bRCheckMB)
1240                 {
1241                     int d;
1242
1243                     for (d = 0; (d < DIM && bUse); d++)
1244                     {
1245                         /* Check if the cg_cm distance falls within
1246                          * the cut-off to avoid possible multiple
1247                          * assignments of bonded interactions.
1248                          */
1249                         if (rcheck[d] && k_plus[d]
1250                             && dd_dist2(pbc_null, cg_cm, tiatoms[k_zero[d]], tiatoms[k_plus[d]]) >= rc2)
1251                         {
1252                             bUse = FALSE;
1253                         }
1254                     }
1255                 }
1256             }
1257             if (bUse)
1258             {
1259                 /* Add this interaction to the local topology */
1260                 idef->il[ftype].push_back(tiatoms[0], nral, tiatoms + 1);
1261                 /* Sum so we can check in global_stat
1262                  * if we have everything.
1263                  */
1264                 if (bBCheck || !(interaction_function[ftype].flags & IF_LIMZERO))
1265                 {
1266                     (*nbonded_local)++;
1267                 }
1268             }
1269         }
1270         j += 1 + nral_rt(ftype);
1271     }
1272 }
1273
1274 /*! \brief This function looks up and assigns bonded interactions for zone iz.
1275  *
1276  * With thread parallelizing each thread acts on a different atom range:
1277  * at_start to at_end.
1278  */
1279 static int make_bondeds_zone(gmx_domdec_t*                      dd,
1280                              const gmx_domdec_zones_t*          zones,
1281                              const std::vector<gmx_molblock_t>& molb,
1282                              gmx_bool                           bRCheckMB,
1283                              ivec                               rcheck,
1284                              gmx_bool                           bRCheck2B,
1285                              real                               rc2,
1286                              t_pbc*                             pbc_null,
1287                              rvec*                              cg_cm,
1288                              const t_iparams*                   ip_in,
1289                              InteractionDefinitions*            idef,
1290                              int                                izone,
1291                              const gmx::Range<int>&             atomRange)
1292 {
1293     int                mb, mt, mol, i_mol;
1294     gmx_bool           bBCheck;
1295     gmx_reverse_top_t* rt;
1296     int                nbonded_local;
1297
1298     rt = dd->reverse_top;
1299
1300     bBCheck = rt->bBCheck;
1301
1302     nbonded_local = 0;
1303
1304     for (int i : atomRange)
1305     {
1306         /* Get the global atom number */
1307         const int i_gl = dd->globalAtomIndices[i];
1308         global_atomnr_to_moltype_ind(rt, i_gl, &mb, &mt, &mol, &i_mol);
1309         /* Check all intramolecular interactions assigned to this atom */
1310         gmx::ArrayRef<const int>     index = rt->ril_mt[mt].index;
1311         gmx::ArrayRef<const t_iatom> rtil  = rt->ril_mt[mt].il;
1312
1313         check_assign_interactions_atom(i,
1314                                        i_gl,
1315                                        mol,
1316                                        i_mol,
1317                                        rt->ril_mt[mt].numAtomsInMolecule,
1318                                        index,
1319                                        rtil,
1320                                        FALSE,
1321                                        index[i_mol],
1322                                        index[i_mol + 1],
1323                                        dd,
1324                                        zones,
1325                                        &molb[mb],
1326                                        bRCheckMB,
1327                                        rcheck,
1328                                        bRCheck2B,
1329                                        rc2,
1330                                        pbc_null,
1331                                        cg_cm,
1332                                        ip_in,
1333                                        idef,
1334                                        izone,
1335                                        bBCheck,
1336                                        &nbonded_local);
1337
1338
1339         if (rt->bIntermolecularInteractions)
1340         {
1341             /* Check all intermolecular interactions assigned to this atom */
1342             index = rt->ril_intermol.index;
1343             rtil  = rt->ril_intermol.il;
1344
1345             check_assign_interactions_atom(i,
1346                                            i_gl,
1347                                            mol,
1348                                            i_mol,
1349                                            rt->ril_mt[mt].numAtomsInMolecule,
1350                                            index,
1351                                            rtil,
1352                                            TRUE,
1353                                            index[i_gl],
1354                                            index[i_gl + 1],
1355                                            dd,
1356                                            zones,
1357                                            &molb[mb],
1358                                            bRCheckMB,
1359                                            rcheck,
1360                                            bRCheck2B,
1361                                            rc2,
1362                                            pbc_null,
1363                                            cg_cm,
1364                                            ip_in,
1365                                            idef,
1366                                            izone,
1367                                            bBCheck,
1368                                            &nbonded_local);
1369         }
1370     }
1371
1372     return nbonded_local;
1373 }
1374
1375 /*! \brief Set the exclusion data for i-zone \p iz */
1376 static void make_exclusions_zone(gmx_domdec_t*                     dd,
1377                                  gmx_domdec_zones_t*               zones,
1378                                  const std::vector<gmx_moltype_t>& moltype,
1379                                  const int*                        cginfo,
1380                                  ListOfLists<int>*                 lexcls,
1381                                  int                               iz,
1382                                  int                               at_start,
1383                                  int                               at_end,
1384                                  const gmx::ArrayRef<const int>    intermolecularExclusionGroup)
1385 {
1386     const gmx_ga2la_t& ga2la = *dd->ga2la;
1387
1388     const auto& jAtomRange = zones->iZones[iz].jAtomRange;
1389
1390     const gmx::index oldNumLists = lexcls->ssize();
1391
1392     std::vector<int> exclusionsForAtom;
1393     for (int at = at_start; at < at_end; at++)
1394     {
1395         exclusionsForAtom.clear();
1396
1397         if (GET_CGINFO_EXCL_INTER(cginfo[at]))
1398         {
1399             int a_gl, mb, mt, mol, a_mol;
1400
1401             /* Copy the exclusions from the global top */
1402             a_gl = dd->globalAtomIndices[at];
1403             global_atomnr_to_moltype_ind(dd->reverse_top, a_gl, &mb, &mt, &mol, &a_mol);
1404             const auto excls = moltype[mt].excls[a_mol];
1405             for (const int aj_mol : excls)
1406             {
1407                 if (const auto* jEntry = ga2la.find(a_gl + aj_mol - a_mol))
1408                 {
1409                     /* This check is not necessary, but it can reduce
1410                      * the number of exclusions in the list, which in turn
1411                      * can speed up the pair list construction a bit.
1412                      */
1413                     if (jAtomRange.isInRange(jEntry->la))
1414                     {
1415                         exclusionsForAtom.push_back(jEntry->la);
1416                     }
1417                 }
1418             }
1419         }
1420
1421         bool isExcludedAtom = !intermolecularExclusionGroup.empty()
1422                               && std::find(intermolecularExclusionGroup.begin(),
1423                                            intermolecularExclusionGroup.end(),
1424                                            dd->globalAtomIndices[at])
1425                                          != intermolecularExclusionGroup.end();
1426
1427         if (isExcludedAtom)
1428         {
1429             for (int qmAtomGlobalIndex : intermolecularExclusionGroup)
1430             {
1431                 if (const auto* entry = dd->ga2la->find(qmAtomGlobalIndex))
1432                 {
1433                     exclusionsForAtom.push_back(entry->la);
1434                 }
1435             }
1436         }
1437
1438         /* Append the exclusions for this atom to the topology */
1439         lexcls->pushBack(exclusionsForAtom);
1440     }
1441
1442     GMX_RELEASE_ASSERT(
1443             lexcls->ssize() - oldNumLists == at_end - at_start,
1444             "The number of exclusion list should match the number of atoms in the range");
1445 }
1446
1447 /*! \brief Generate and store all required local bonded interactions in \p idef and local exclusions in \p lexcls */
1448 static int make_local_bondeds_excls(gmx_domdec_t*           dd,
1449                                     gmx_domdec_zones_t*     zones,
1450                                     const gmx_mtop_t*       mtop,
1451                                     const int*              cginfo,
1452                                     gmx_bool                bRCheckMB,
1453                                     ivec                    rcheck,
1454                                     gmx_bool                bRCheck2B,
1455                                     real                    rc,
1456                                     t_pbc*                  pbc_null,
1457                                     rvec*                   cg_cm,
1458                                     InteractionDefinitions* idef,
1459                                     ListOfLists<int>*       lexcls,
1460                                     int*                    excl_count)
1461 {
1462     int                nzone_bondeds;
1463     int                cg0, cg1;
1464     real               rc2;
1465     int                nbonded_local;
1466     gmx_reverse_top_t* rt;
1467
1468     if (dd->reverse_top->bInterAtomicInteractions)
1469     {
1470         nzone_bondeds = zones->n;
1471     }
1472     else
1473     {
1474         /* Only single charge group (or atom) molecules, so interactions don't
1475          * cross zone boundaries and we only need to assign in the home zone.
1476          */
1477         nzone_bondeds = 1;
1478     }
1479
1480     /* We only use exclusions from i-zones to i- and j-zones */
1481     const int numIZonesForExclusions = (dd->haveExclusions ? zones->iZones.size() : 0);
1482
1483     rt = dd->reverse_top;
1484
1485     rc2 = rc * rc;
1486
1487     /* Clear the counts */
1488     idef->clear();
1489     nbonded_local = 0;
1490
1491     lexcls->clear();
1492     *excl_count = 0;
1493
1494     for (int izone = 0; izone < nzone_bondeds; izone++)
1495     {
1496         cg0 = zones->cg_range[izone];
1497         cg1 = zones->cg_range[izone + 1];
1498
1499         const int numThreads = rt->th_work.size();
1500 #pragma omp parallel for num_threads(numThreads) schedule(static)
1501         for (int thread = 0; thread < numThreads; thread++)
1502         {
1503             try
1504             {
1505                 int                     cg0t, cg1t;
1506                 InteractionDefinitions* idef_t;
1507
1508                 cg0t = cg0 + ((cg1 - cg0) * thread) / numThreads;
1509                 cg1t = cg0 + ((cg1 - cg0) * (thread + 1)) / numThreads;
1510
1511                 if (thread == 0)
1512                 {
1513                     idef_t = idef;
1514                 }
1515                 else
1516                 {
1517                     idef_t = &rt->th_work[thread].idef;
1518                     idef_t->clear();
1519                 }
1520
1521                 rt->th_work[thread].nbonded = make_bondeds_zone(dd,
1522                                                                 zones,
1523                                                                 mtop->molblock,
1524                                                                 bRCheckMB,
1525                                                                 rcheck,
1526                                                                 bRCheck2B,
1527                                                                 rc2,
1528                                                                 pbc_null,
1529                                                                 cg_cm,
1530                                                                 idef->iparams.data(),
1531                                                                 idef_t,
1532                                                                 izone,
1533                                                                 gmx::Range<int>(cg0t, cg1t));
1534
1535                 if (izone < numIZonesForExclusions)
1536                 {
1537                     ListOfLists<int>* excl_t;
1538                     if (thread == 0)
1539                     {
1540                         // Thread 0 stores exclusions directly in the final storage
1541                         excl_t = lexcls;
1542                     }
1543                     else
1544                     {
1545                         // Threads > 0 store in temporary storage, starting at list index 0
1546                         excl_t = &rt->th_work[thread].excl;
1547                         excl_t->clear();
1548                     }
1549
1550                     /* No charge groups and no distance check required */
1551                     make_exclusions_zone(
1552                             dd, zones, mtop->moltype, cginfo, excl_t, izone, cg0t, cg1t, mtop->intermolecularExclusionGroup);
1553                 }
1554             }
1555             GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
1556         }
1557
1558         if (rt->th_work.size() > 1)
1559         {
1560             combine_idef(idef, rt->th_work);
1561         }
1562
1563         for (const thread_work_t& th_work : rt->th_work)
1564         {
1565             nbonded_local += th_work.nbonded;
1566         }
1567
1568         if (izone < numIZonesForExclusions)
1569         {
1570             for (std::size_t th = 1; th < rt->th_work.size(); th++)
1571             {
1572                 lexcls->appendListOfLists(rt->th_work[th].excl);
1573             }
1574             for (const thread_work_t& th_work : rt->th_work)
1575             {
1576                 *excl_count += th_work.excl_count;
1577             }
1578         }
1579     }
1580
1581     if (debug)
1582     {
1583         fprintf(debug, "We have %d exclusions, check count %d\n", lexcls->numElements(), *excl_count);
1584     }
1585
1586     return nbonded_local;
1587 }
1588
1589 void dd_make_local_top(gmx_domdec_t*       dd,
1590                        gmx_domdec_zones_t* zones,
1591                        int                 npbcdim,
1592                        matrix              box,
1593                        rvec                cellsize_min,
1594                        const ivec          npulse,
1595                        t_forcerec*         fr,
1596                        rvec*               cgcm_or_x,
1597                        const gmx_mtop_t&   mtop,
1598                        gmx_localtop_t*     ltop)
1599 {
1600     gmx_bool bRCheckMB, bRCheck2B;
1601     real     rc = -1;
1602     ivec     rcheck;
1603     int      d, nexcl;
1604     t_pbc    pbc, *pbc_null = nullptr;
1605
1606     if (debug)
1607     {
1608         fprintf(debug, "Making local topology\n");
1609     }
1610
1611     bRCheckMB = FALSE;
1612     bRCheck2B = FALSE;
1613
1614     if (dd->reverse_top->bInterAtomicInteractions)
1615     {
1616         /* We need to check to which cell bondeds should be assigned */
1617         rc = dd_cutoff_twobody(dd);
1618         if (debug)
1619         {
1620             fprintf(debug, "Two-body bonded cut-off distance is %g\n", rc);
1621         }
1622
1623         /* Should we check cg_cm distances when assigning bonded interactions? */
1624         for (d = 0; d < DIM; d++)
1625         {
1626             rcheck[d] = FALSE;
1627             /* Only need to check for dimensions where the part of the box
1628              * that is not communicated is smaller than the cut-off.
1629              */
1630             if (d < npbcdim && dd->numCells[d] > 1
1631                 && (dd->numCells[d] - npulse[d]) * cellsize_min[d] < 2 * rc)
1632             {
1633                 if (dd->numCells[d] == 2)
1634                 {
1635                     rcheck[d] = TRUE;
1636                     bRCheckMB = TRUE;
1637                 }
1638                 /* Check for interactions between two atoms,
1639                  * where we can allow interactions up to the cut-off,
1640                  * instead of up to the smallest cell dimension.
1641                  */
1642                 bRCheck2B = TRUE;
1643             }
1644             if (debug)
1645             {
1646                 fprintf(debug,
1647                         "dim %d cellmin %f bonded rcheck[%d] = %d, bRCheck2B = %s\n",
1648                         d,
1649                         cellsize_min[d],
1650                         d,
1651                         rcheck[d],
1652                         gmx::boolToString(bRCheck2B));
1653             }
1654         }
1655         if (bRCheckMB || bRCheck2B)
1656         {
1657             if (fr->bMolPBC)
1658             {
1659                 pbc_null = set_pbc_dd(&pbc, fr->pbcType, dd->numCells, TRUE, box);
1660             }
1661             else
1662             {
1663                 pbc_null = nullptr;
1664             }
1665         }
1666     }
1667
1668     dd->nbonded_local = make_local_bondeds_excls(dd,
1669                                                  zones,
1670                                                  &mtop,
1671                                                  fr->cginfo.data(),
1672                                                  bRCheckMB,
1673                                                  rcheck,
1674                                                  bRCheck2B,
1675                                                  rc,
1676                                                  pbc_null,
1677                                                  cgcm_or_x,
1678                                                  &ltop->idef,
1679                                                  &ltop->excls,
1680                                                  &nexcl);
1681
1682     /* The ilist is not sorted yet,
1683      * we can only do this when we have the charge arrays.
1684      */
1685     ltop->idef.ilsort = ilsortUNKNOWN;
1686 }
1687
1688 void dd_sort_local_top(gmx_domdec_t* dd, const t_mdatoms* mdatoms, gmx_localtop_t* ltop)
1689 {
1690     if (dd->reverse_top->ilsort == ilsortNO_FE)
1691     {
1692         ltop->idef.ilsort = ilsortNO_FE;
1693     }
1694     else
1695     {
1696         gmx_sort_ilist_fe(&ltop->idef, mdatoms->chargeA, mdatoms->chargeB);
1697     }
1698 }
1699
1700 void dd_init_local_state(gmx_domdec_t* dd, const t_state* state_global, t_state* state_local)
1701 {
1702     int buf[NITEM_DD_INIT_LOCAL_STATE];
1703
1704     if (DDMASTER(dd))
1705     {
1706         buf[0] = state_global->flags;
1707         buf[1] = state_global->ngtc;
1708         buf[2] = state_global->nnhpres;
1709         buf[3] = state_global->nhchainlength;
1710         buf[4] = state_global->dfhist ? state_global->dfhist->nlambda : 0;
1711     }
1712     dd_bcast(dd, NITEM_DD_INIT_LOCAL_STATE * sizeof(int), buf);
1713
1714     init_gtc_state(state_local, buf[1], buf[2], buf[3]);
1715     init_dfhist_state(state_local, buf[4]);
1716     state_local->flags = buf[0];
1717 }
1718
1719 /*! \brief Check if a link is stored in \p link between charge groups \p cg_gl and \p cg_gl_j and if not so, store a link */
1720 static void check_link(t_blocka* link, int cg_gl, int cg_gl_j)
1721 {
1722     int      k;
1723     gmx_bool bFound;
1724
1725     bFound = FALSE;
1726     for (k = link->index[cg_gl]; k < link->index[cg_gl + 1]; k++)
1727     {
1728         GMX_RELEASE_ASSERT(link->a, "Inconsistent NULL pointer while making charge-group links");
1729         if (link->a[k] == cg_gl_j)
1730         {
1731             bFound = TRUE;
1732         }
1733     }
1734     if (!bFound)
1735     {
1736         GMX_RELEASE_ASSERT(link->a || link->index[cg_gl + 1] + 1 > link->nalloc_a,
1737                            "Inconsistent allocation of link");
1738         /* Add this charge group link */
1739         if (link->index[cg_gl + 1] + 1 > link->nalloc_a)
1740         {
1741             link->nalloc_a = over_alloc_large(link->index[cg_gl + 1] + 1);
1742             srenew(link->a, link->nalloc_a);
1743         }
1744         link->a[link->index[cg_gl + 1]] = cg_gl_j;
1745         link->index[cg_gl + 1]++;
1746     }
1747 }
1748
1749 t_blocka* makeBondedLinks(const gmx_mtop_t& mtop, gmx::ArrayRef<cginfo_mb_t> cginfo_mb)
1750 {
1751     t_blocka*    link;
1752     cginfo_mb_t* cgi_mb;
1753
1754     /* For each atom make a list of other atoms in the system
1755      * that a linked to it via bonded interactions
1756      * which are also stored in reverse_top.
1757      */
1758
1759     reverse_ilist_t ril_intermol;
1760     if (mtop.bIntermolecularInteractions)
1761     {
1762         t_atoms atoms;
1763
1764         atoms.nr   = mtop.natoms;
1765         atoms.atom = nullptr;
1766
1767         GMX_RELEASE_ASSERT(mtop.intermolecular_ilist,
1768                            "We should have an ilist when intermolecular interactions are on");
1769
1770         make_reverse_ilist(*mtop.intermolecular_ilist, &atoms, FALSE, FALSE, FALSE, TRUE, &ril_intermol);
1771     }
1772
1773     snew(link, 1);
1774     snew(link->index, mtop.natoms + 1);
1775     link->nalloc_a = 0;
1776     link->a        = nullptr;
1777
1778     link->index[0] = 0;
1779     int cg_offset  = 0;
1780     int ncgi       = 0;
1781     for (size_t mb = 0; mb < mtop.molblock.size(); mb++)
1782     {
1783         const gmx_molblock_t& molb = mtop.molblock[mb];
1784         if (molb.nmol == 0)
1785         {
1786             continue;
1787         }
1788         const gmx_moltype_t& molt = mtop.moltype[molb.type];
1789         /* Make a reverse ilist in which the interactions are linked
1790          * to all atoms, not only the first atom as in gmx_reverse_top.
1791          * The constraints are discarded here.
1792          */
1793         reverse_ilist_t ril;
1794         make_reverse_ilist(molt.ilist, &molt.atoms, FALSE, FALSE, FALSE, TRUE, &ril);
1795
1796         cgi_mb = &cginfo_mb[mb];
1797
1798         int mol;
1799         for (mol = 0; mol < (mtop.bIntermolecularInteractions ? molb.nmol : 1); mol++)
1800         {
1801             for (int a = 0; a < molt.atoms.nr; a++)
1802             {
1803                 int cg_gl              = cg_offset + a;
1804                 link->index[cg_gl + 1] = link->index[cg_gl];
1805                 int i                  = ril.index[a];
1806                 while (i < ril.index[a + 1])
1807                 {
1808                     int ftype = ril.il[i++];
1809                     int nral  = NRAL(ftype);
1810                     /* Skip the ifunc index */
1811                     i++;
1812                     for (int j = 0; j < nral; j++)
1813                     {
1814                         int aj = ril.il[i + j];
1815                         if (aj != a)
1816                         {
1817                             check_link(link, cg_gl, cg_offset + aj);
1818                         }
1819                     }
1820                     i += nral_rt(ftype);
1821                 }
1822
1823                 if (mtop.bIntermolecularInteractions)
1824                 {
1825                     int i = ril_intermol.index[cg_gl];
1826                     while (i < ril_intermol.index[cg_gl + 1])
1827                     {
1828                         int ftype = ril_intermol.il[i++];
1829                         int nral  = NRAL(ftype);
1830                         /* Skip the ifunc index */
1831                         i++;
1832                         for (int j = 0; j < nral; j++)
1833                         {
1834                             /* Here we assume we have no charge groups;
1835                              * this has been checked above.
1836                              */
1837                             int aj = ril_intermol.il[i + j];
1838                             check_link(link, cg_gl, aj);
1839                         }
1840                         i += nral_rt(ftype);
1841                     }
1842                 }
1843                 if (link->index[cg_gl + 1] - link->index[cg_gl] > 0)
1844                 {
1845                     SET_CGINFO_BOND_INTER(cgi_mb->cginfo[a]);
1846                     ncgi++;
1847                 }
1848             }
1849
1850             cg_offset += molt.atoms.nr;
1851         }
1852         int nlink_mol = link->index[cg_offset] - link->index[cg_offset - molt.atoms.nr];
1853
1854         if (debug)
1855         {
1856             fprintf(debug,
1857                     "molecule type '%s' %d atoms has %d atom links through bonded interac.\n",
1858                     *molt.name,
1859                     molt.atoms.nr,
1860                     nlink_mol);
1861         }
1862
1863         if (molb.nmol > mol)
1864         {
1865             /* Copy the data for the rest of the molecules in this block */
1866             link->nalloc_a += (molb.nmol - mol) * nlink_mol;
1867             srenew(link->a, link->nalloc_a);
1868             for (; mol < molb.nmol; mol++)
1869             {
1870                 for (int a = 0; a < molt.atoms.nr; a++)
1871                 {
1872                     int cg_gl              = cg_offset + a;
1873                     link->index[cg_gl + 1] = link->index[cg_gl + 1 - molt.atoms.nr] + nlink_mol;
1874                     for (int j = link->index[cg_gl]; j < link->index[cg_gl + 1]; j++)
1875                     {
1876                         link->a[j] = link->a[j - nlink_mol] + molt.atoms.nr;
1877                     }
1878                     if (link->index[cg_gl + 1] - link->index[cg_gl] > 0
1879                         && cg_gl - cgi_mb->cg_start < cgi_mb->cg_mod)
1880                     {
1881                         SET_CGINFO_BOND_INTER(cgi_mb->cginfo[cg_gl - cgi_mb->cg_start]);
1882                         ncgi++;
1883                     }
1884                 }
1885                 cg_offset += molt.atoms.nr;
1886             }
1887         }
1888     }
1889
1890     if (debug)
1891     {
1892         fprintf(debug, "Of the %d atoms %d are linked via bonded interactions\n", mtop.natoms, ncgi);
1893     }
1894
1895     return link;
1896 }
1897
1898 typedef struct
1899 {
1900     real r2;
1901     int  ftype;
1902     int  a1;
1903     int  a2;
1904 } bonded_distance_t;
1905
1906 /*! \brief Compare distance^2 \p r2 against the distance in \p bd and if larger store it along with \p ftype and atom indices \p a1 and \p a2 */
1907 static void update_max_bonded_distance(real r2, int ftype, int a1, int a2, bonded_distance_t* bd)
1908 {
1909     if (r2 > bd->r2)
1910     {
1911         bd->r2    = r2;
1912         bd->ftype = ftype;
1913         bd->a1    = a1;
1914         bd->a2    = a2;
1915     }
1916 }
1917
1918 /*! \brief Set the distance, function type and atom indices for the longest distance between atoms of molecule type \p molt for two-body and multi-body bonded interactions */
1919 static void bonded_cg_distance_mol(const gmx_moltype_t* molt,
1920                                    gmx_bool             bBCheck,
1921                                    gmx_bool             bExcl,
1922                                    ArrayRef<const RVec> x,
1923                                    bonded_distance_t*   bd_2b,
1924                                    bonded_distance_t*   bd_mb)
1925 {
1926     for (int ftype = 0; ftype < F_NRE; ftype++)
1927     {
1928         if (dd_check_ftype(ftype, bBCheck, FALSE, FALSE))
1929         {
1930             const auto& il   = molt->ilist[ftype];
1931             int         nral = NRAL(ftype);
1932             if (nral > 1)
1933             {
1934                 for (int i = 0; i < il.size(); i += 1 + nral)
1935                 {
1936                     for (int ai = 0; ai < nral; ai++)
1937                     {
1938                         int atomI = il.iatoms[i + 1 + ai];
1939                         for (int aj = ai + 1; aj < nral; aj++)
1940                         {
1941                             int atomJ = il.iatoms[i + 1 + aj];
1942                             if (atomI != atomJ)
1943                             {
1944                                 real rij2 = distance2(x[atomI], x[atomJ]);
1945
1946                                 update_max_bonded_distance(
1947                                         rij2, ftype, atomI, atomJ, (nral == 2) ? bd_2b : bd_mb);
1948                             }
1949                         }
1950                     }
1951                 }
1952             }
1953         }
1954     }
1955     if (bExcl)
1956     {
1957         const auto& excls = molt->excls;
1958         for (gmx::index ai = 0; ai < excls.ssize(); ai++)
1959         {
1960             for (const int aj : excls[ai])
1961             {
1962                 if (ai != aj)
1963                 {
1964                     real rij2 = distance2(x[ai], x[aj]);
1965
1966                     /* There is no function type for exclusions, use -1 */
1967                     update_max_bonded_distance(rij2, -1, ai, aj, bd_2b);
1968                 }
1969             }
1970         }
1971     }
1972 }
1973
1974 /*! \brief Set the distance, function type and atom indices for the longest atom distance involved in intermolecular interactions for two-body and multi-body bonded interactions */
1975 static void bonded_distance_intermol(const InteractionLists& ilists_intermol,
1976                                      gmx_bool                bBCheck,
1977                                      ArrayRef<const RVec>    x,
1978                                      PbcType                 pbcType,
1979                                      const matrix            box,
1980                                      bonded_distance_t*      bd_2b,
1981                                      bonded_distance_t*      bd_mb)
1982 {
1983     t_pbc pbc;
1984
1985     set_pbc(&pbc, pbcType, box);
1986
1987     for (int ftype = 0; ftype < F_NRE; ftype++)
1988     {
1989         if (dd_check_ftype(ftype, bBCheck, FALSE, FALSE))
1990         {
1991             const auto& il   = ilists_intermol[ftype];
1992             int         nral = NRAL(ftype);
1993
1994             /* No nral>1 check here, since intermol interactions always
1995              * have nral>=2 (and the code is also correct for nral=1).
1996              */
1997             for (int i = 0; i < il.size(); i += 1 + nral)
1998             {
1999                 for (int ai = 0; ai < nral; ai++)
2000                 {
2001                     int atom_i = il.iatoms[i + 1 + ai];
2002
2003                     for (int aj = ai + 1; aj < nral; aj++)
2004                     {
2005                         rvec dx;
2006                         real rij2;
2007
2008                         int atom_j = il.iatoms[i + 1 + aj];
2009
2010                         pbc_dx(&pbc, x[atom_i], x[atom_j], dx);
2011
2012                         rij2 = norm2(dx);
2013
2014                         update_max_bonded_distance(rij2, ftype, atom_i, atom_j, (nral == 2) ? bd_2b : bd_mb);
2015                     }
2016                 }
2017             }
2018         }
2019     }
2020 }
2021
2022 //! Returns whether \p molt has at least one virtual site
2023 static bool moltypeHasVsite(const gmx_moltype_t& molt)
2024 {
2025     bool hasVsite = false;
2026     for (int i = 0; i < F_NRE; i++)
2027     {
2028         if ((interaction_function[i].flags & IF_VSITE) && !molt.ilist[i].empty())
2029         {
2030             hasVsite = true;
2031         }
2032     }
2033
2034     return hasVsite;
2035 }
2036
2037 //! Returns coordinates not broken over PBC for a molecule
2038 static void getWholeMoleculeCoordinates(const gmx_moltype_t*  molt,
2039                                         const gmx_ffparams_t* ffparams,
2040                                         PbcType               pbcType,
2041                                         t_graph*              graph,
2042                                         const matrix          box,
2043                                         ArrayRef<const RVec>  x,
2044                                         ArrayRef<RVec>        xs)
2045 {
2046     int n, i;
2047
2048     if (pbcType != PbcType::No)
2049     {
2050         mk_mshift(nullptr, graph, pbcType, box, as_rvec_array(x.data()));
2051
2052         shift_x(graph, box, as_rvec_array(x.data()), as_rvec_array(xs.data()));
2053         /* By doing an extra mk_mshift the molecules that are broken
2054          * because they were e.g. imported from another software
2055          * will be made whole again. Such are the healing powers
2056          * of GROMACS.
2057          */
2058         mk_mshift(nullptr, graph, pbcType, box, as_rvec_array(xs.data()));
2059     }
2060     else
2061     {
2062         /* We copy the coordinates so the original coordinates remain
2063          * unchanged, just to be 100% sure that we do not affect
2064          * binary reproducibility of simulations.
2065          */
2066         n = molt->atoms.nr;
2067         for (i = 0; i < n; i++)
2068         {
2069             copy_rvec(x[i], xs[i]);
2070         }
2071     }
2072
2073     if (moltypeHasVsite(*molt))
2074     {
2075         gmx::constructVirtualSites(xs, ffparams->iparams, molt->ilist);
2076     }
2077 }
2078
2079 void dd_bonded_cg_distance(const gmx::MDLogger& mdlog,
2080                            const gmx_mtop_t*    mtop,
2081                            const t_inputrec*    ir,
2082                            ArrayRef<const RVec> x,
2083                            const matrix         box,
2084                            gmx_bool             bBCheck,
2085                            real*                r_2b,
2086                            real*                r_mb)
2087 {
2088     gmx_bool          bExclRequired;
2089     int               at_offset;
2090     bonded_distance_t bd_2b = { 0, -1, -1, -1 };
2091     bonded_distance_t bd_mb = { 0, -1, -1, -1 };
2092
2093     bExclRequired = inputrecExclForces(ir);
2094
2095     *r_2b     = 0;
2096     *r_mb     = 0;
2097     at_offset = 0;
2098     for (const gmx_molblock_t& molb : mtop->molblock)
2099     {
2100         const gmx_moltype_t& molt = mtop->moltype[molb.type];
2101         if (molt.atoms.nr == 1 || molb.nmol == 0)
2102         {
2103             at_offset += molb.nmol * molt.atoms.nr;
2104         }
2105         else
2106         {
2107             t_graph graph;
2108             if (ir->pbcType != PbcType::No)
2109             {
2110                 graph = mk_graph_moltype(molt);
2111             }
2112
2113             std::vector<RVec> xs(molt.atoms.nr);
2114             for (int mol = 0; mol < molb.nmol; mol++)
2115             {
2116                 getWholeMoleculeCoordinates(&molt,
2117                                             &mtop->ffparams,
2118                                             ir->pbcType,
2119                                             &graph,
2120                                             box,
2121                                             x.subArray(at_offset, molt.atoms.nr),
2122                                             xs);
2123
2124                 bonded_distance_t bd_mol_2b = { 0, -1, -1, -1 };
2125                 bonded_distance_t bd_mol_mb = { 0, -1, -1, -1 };
2126
2127                 bonded_cg_distance_mol(&molt, bBCheck, bExclRequired, xs, &bd_mol_2b, &bd_mol_mb);
2128
2129                 /* Process the mol data adding the atom index offset */
2130                 update_max_bonded_distance(bd_mol_2b.r2,
2131                                            bd_mol_2b.ftype,
2132                                            at_offset + bd_mol_2b.a1,
2133                                            at_offset + bd_mol_2b.a2,
2134                                            &bd_2b);
2135                 update_max_bonded_distance(bd_mol_mb.r2,
2136                                            bd_mol_mb.ftype,
2137                                            at_offset + bd_mol_mb.a1,
2138                                            at_offset + bd_mol_mb.a2,
2139                                            &bd_mb);
2140
2141                 at_offset += molt.atoms.nr;
2142             }
2143         }
2144     }
2145
2146     if (mtop->bIntermolecularInteractions)
2147     {
2148         GMX_RELEASE_ASSERT(mtop->intermolecular_ilist,
2149                            "We should have an ilist when intermolecular interactions are on");
2150
2151         bonded_distance_intermol(*mtop->intermolecular_ilist, bBCheck, x, ir->pbcType, box, &bd_2b, &bd_mb);
2152     }
2153
2154     *r_2b = sqrt(bd_2b.r2);
2155     *r_mb = sqrt(bd_mb.r2);
2156
2157     if (*r_2b > 0 || *r_mb > 0)
2158     {
2159         GMX_LOG(mdlog.info).appendText("Initial maximum distances in bonded interactions:");
2160         if (*r_2b > 0)
2161         {
2162             GMX_LOG(mdlog.info)
2163                     .appendTextFormatted(
2164                             "    two-body bonded interactions: %5.3f nm, %s, atoms %d %d",
2165                             *r_2b,
2166                             (bd_2b.ftype >= 0) ? interaction_function[bd_2b.ftype].longname : "Exclusion",
2167                             bd_2b.a1 + 1,
2168                             bd_2b.a2 + 1);
2169         }
2170         if (*r_mb > 0)
2171         {
2172             GMX_LOG(mdlog.info)
2173                     .appendTextFormatted(
2174                             "  multi-body bonded interactions: %5.3f nm, %s, atoms %d %d",
2175                             *r_mb,
2176                             interaction_function[bd_mb.ftype].longname,
2177                             bd_mb.a1 + 1,
2178                             bd_mb.a2 + 1);
2179         }
2180     }
2181 }