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