src/gromacs/domdec/domdec_topology.cpp

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