src/gromacs/domdec/domdec_topology.cpp

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