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