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