src/gromacs/domdec/domdec_topology.cpp

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