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