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