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