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