src/gromacs/domdec/domdec_topology.cpp

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