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 <optional>
  56 #include <string>
  57
  58 #include "gromacs/domdec/domdec.h"
  59 #include "gromacs/domdec/domdec_network.h"
  60 #include "gromacs/domdec/ga2la.h"
  61 #include "gromacs/domdec/localtopologychecker.h"
  62 #include "gromacs/domdec/options.h"
  63 #include "gromacs/domdec/reversetopology.h"
  64 #include "gromacs/gmxlib/network.h"
  65 #include "gromacs/math/vec.h"
  66 #include "gromacs/mdlib/forcerec.h"
  67 #include "gromacs/mdlib/gmx_omp_nthreads.h"
  68 #include "gromacs/mdlib/vsite.h"
  69 #include "gromacs/mdtypes/commrec.h"
  70 #include "gromacs/mdtypes/forcerec.h"
  71 #include "gromacs/mdtypes/inputrec.h"
  72 #include "gromacs/mdtypes/md_enums.h"
  73 #include "gromacs/mdtypes/mdatom.h"
  74 #include "gromacs/mdtypes/state.h"
  75 #include "gromacs/pbcutil/mshift.h"
  76 #include "gromacs/pbcutil/pbc.h"
  77 #include "gromacs/topology/mtop_util.h"
  78 #include "gromacs/topology/topsort.h"
  79 #include "gromacs/utility/cstringutil.h"
  80 #include "gromacs/utility/exceptions.h"
  81 #include "gromacs/utility/fatalerror.h"
  82 #include "gromacs/utility/gmxassert.h"
  83 #include "gromacs/utility/logger.h"
  84 #include "gromacs/utility/smalloc.h"
  85 #include "gromacs/utility/strconvert.h"
  86 #include "gromacs/utility/stringstream.h"
  87 #include "gromacs/utility/stringutil.h"
  88 #include "gromacs/utility/textwriter.h"
  89
  90 #include "domdec_constraints.h"
  91 #include "domdec_internal.h"
  92 #include "domdec_vsite.h"
  93 #include "dump.h"
  94
  95 using gmx::ArrayRef;
  96 using gmx::DDBondedChecking;
  97 using gmx::ListOfLists;
  98 using gmx::RVec;
  99
 100 /*! \brief The number of integer item in the local state, used for broadcasting of the state */
 101 #define NITEM_DD_INIT_LOCAL_STATE 5
 102
 103 /*! \brief Store a vsite interaction at the end of \p il
 104  *
 105  * This routine is very similar to add_ifunc, but vsites interactions
 106  * have more work to do than other kinds of interactions, and the
 107  * complex way nral (and thus vector contents) depends on ftype
 108  * confuses static analysis tools unless we fuse the vsite
 109  * atom-indexing organization code with the ifunc-adding code, so that
 110  * they can see that nral is the same value. */
 111 static inline void add_ifunc_for_vsites(t_iatom*           tiatoms,
 112                                         const gmx_ga2la_t& ga2la,
 113                                         int                nral,
 114                                         gmx_bool           bHomeA,
 115                                         int                a,
 116                                         int                a_gl,
 117                                         int                a_mol,
 118                                         const t_iatom*     iatoms,
 119                                         InteractionList*   il)
 120 {
 121     /* Copy the type */
 122     tiatoms[0] = iatoms[0];
 123
 124     if (bHomeA)
 125     {
 126         /* We know the local index of the first atom */
 127         tiatoms[1] = a;
 128     }
 129     else
 130     {
 131         /* Convert later in make_local_vsites */
 132         tiatoms[1] = -a_gl - 1;
 133     }
 134
 135     GMX_ASSERT(nral >= 2 && nral <= 5, "Invalid nral for vsites");
 136     for (int k = 2; k < 1 + nral; k++)
 137     {
 138         int ak_gl = a_gl + iatoms[k] - a_mol;
 139         if (const int* homeIndex = ga2la.findHome(ak_gl))
 140         {
 141             tiatoms[k] = *homeIndex;
 142         }
 143         else
 144         {
 145             /* Copy the global index, convert later in make_local_vsites */
 146             tiatoms[k] = -(ak_gl + 1);
 147         }
 148         // Note that ga2la_get_home always sets the third parameter if
 149         // it returns TRUE
 150     }
 151     il->push_back(tiatoms[0], nral, tiatoms + 1);
 152 }
 153
 154 /*! \brief Store a position restraint in idef and iatoms, complex because the parameters are different for each entry */
 155 static void add_posres(int                     mol,
 156                        int                     a_mol,
 157                        int                     numAtomsInMolecule,
 158                        const gmx_molblock_t&   molb,
 159                        gmx::ArrayRef<int>      iatoms,
 160                        const t_iparams*        ip_in,
 161                        InteractionDefinitions* idef)
 162 {
 163     /* This position restraint has not been added yet,
 164      * so it's index is the current number of position restraints.
 165      */
 166     const int n = idef->il[F_POSRES].size() / 2;
 167
 168     /* Get the position restraint coordinates from the molblock */
 169     const int a_molb = mol * numAtomsInMolecule + a_mol;
 170     GMX_ASSERT(a_molb < ssize(molb.posres_xA),
 171                "We need a sufficient number of position restraint coordinates");
 172     /* Copy the force constants */
 173     t_iparams ip        = ip_in[iatoms[0]];
 174     ip.posres.pos0A[XX] = molb.posres_xA[a_molb][XX];
 175     ip.posres.pos0A[YY] = molb.posres_xA[a_molb][YY];
 176     ip.posres.pos0A[ZZ] = molb.posres_xA[a_molb][ZZ];
 177     if (!molb.posres_xB.empty())
 178     {
 179         ip.posres.pos0B[XX] = molb.posres_xB[a_molb][XX];
 180         ip.posres.pos0B[YY] = molb.posres_xB[a_molb][YY];
 181         ip.posres.pos0B[ZZ] = molb.posres_xB[a_molb][ZZ];
 182     }
 183     else
 184     {
 185         ip.posres.pos0B[XX] = ip.posres.pos0A[XX];
 186         ip.posres.pos0B[YY] = ip.posres.pos0A[YY];
 187         ip.posres.pos0B[ZZ] = ip.posres.pos0A[ZZ];
 188     }
 189     /* Set the parameter index for idef->iparams_posres */
 190     iatoms[0] = n;
 191     idef->iparams_posres.push_back(ip);
 192
 193     GMX_ASSERT(int(idef->iparams_posres.size()) == n + 1,
 194                "The index of the parameter type should match n");
 195 }
 196
 197 /*! \brief Store a flat-bottomed position restraint in idef and iatoms, complex because the parameters are different for each entry */
 198 static void add_fbposres(int                     mol,
 199                          int                     a_mol,
 200                          int                     numAtomsInMolecule,
 201                          const gmx_molblock_t&   molb,
 202                          gmx::ArrayRef<int>      iatoms,
 203                          const t_iparams*        ip_in,
 204                          InteractionDefinitions* idef)
 205 {
 206     /* This flat-bottom position restraint has not been added yet,
 207      * so it's index is the current number of position restraints.
 208      */
 209     const int n = idef->il[F_FBPOSRES].size() / 2;
 210
 211     /* Get the position restraint coordinats from the molblock */
 212     const int a_molb = mol * numAtomsInMolecule + a_mol;
 213     GMX_ASSERT(a_molb < ssize(molb.posres_xA),
 214                "We need a sufficient number of position restraint coordinates");
 215     /* Copy the force constants */
 216     t_iparams ip = ip_in[iatoms[0]];
 217     /* Take reference positions from A position of normal posres */
 218     ip.fbposres.pos0[XX] = molb.posres_xA[a_molb][XX];
 219     ip.fbposres.pos0[YY] = molb.posres_xA[a_molb][YY];
 220     ip.fbposres.pos0[ZZ] = molb.posres_xA[a_molb][ZZ];
 221
 222     /* Note: no B-type for flat-bottom posres */
 223
 224     /* Set the parameter index for idef->iparams_fbposres */
 225     iatoms[0] = n;
 226     idef->iparams_fbposres.push_back(ip);
 227
 228     GMX_ASSERT(int(idef->iparams_fbposres.size()) == n + 1,
 229                "The index of the parameter type should match n");
 230 }
 231
 232 /*! \brief Store a virtual site interaction, complex because of PBC and recursion */
 233 static void add_vsite(const gmx_ga2la_t&       ga2la,
 234                       gmx::ArrayRef<const int> index,
 235                       gmx::ArrayRef<const int> rtil,
 236                       int                      ftype,
 237                       int                      nral,
 238                       gmx_bool                 bHomeA,
 239                       int                      a,
 240                       int                      a_gl,
 241                       int                      a_mol,
 242                       const t_iatom*           iatoms,
 243                       InteractionDefinitions*  idef)
 244 {
 245     t_iatom tiatoms[1 + MAXATOMLIST];
 246
 247     /* Add this interaction to the local topology */
 248     add_ifunc_for_vsites(tiatoms, ga2la, nral, bHomeA, a, a_gl, a_mol, iatoms, &idef->il[ftype]);
 249
 250     if (iatoms[1 + nral])
 251     {
 252         /* Check for recursion */
 253         for (int k = 2; k < 1 + nral; k++)
 254         {
 255             if ((iatoms[1 + nral] & (2 << k)) && (tiatoms[k] < 0))
 256             {
 257                 /* This construction atoms is a vsite and not a home atom */
 258                 if (gmx_debug_at)
 259                 {
 260                     fprintf(debug,
 261                             "Constructing atom %d of vsite atom %d is a vsite and non-home\n",
 262                             iatoms[k] + 1,
 263                             a_mol + 1);
 264                 }
 265                 /* Find the vsite construction */
 266
 267                 /* Check all interactions assigned to this atom */
 268                 int j = index[iatoms[k]];
 269                 while (j < index[iatoms[k] + 1])
 270                 {
 271                     int ftype_r = rtil[j++];
 272                     int nral_r  = NRAL(ftype_r);
 273                     if (interaction_function[ftype_r].flags & IF_VSITE)
 274                     {
 275                         /* Add this vsite (recursion) */
 276                         add_vsite(ga2la,
 277                                   index,
 278                                   rtil,
 279                                   ftype_r,
 280                                   nral_r,
 281                                   FALSE,
 282                                   -1,
 283                                   a_gl + iatoms[k] - iatoms[1],
 284                                   iatoms[k],
 285                                   rtil.data() + j,
 286                                   idef);
 287                     }
 288                     j += 1 + nral_rt(ftype_r);
 289                 }
 290             }
 291         }
 292     }
 293 }
 294
 295 /*! \brief Returns the squared distance between atoms \p i and \p j */
 296 static real dd_dist2(const t_pbc* pbc_null, ArrayRef<const RVec> coordinates, const int i, const int j)
 297 {
 298     rvec dx;
 299
 300     if (pbc_null)
 301     {
 302         pbc_dx_aiuc(pbc_null, coordinates[i], coordinates[j], dx);
 303     }
 304     else
 305     {
 306         rvec_sub(coordinates[i], coordinates[j], dx);
 307     }
 308
 309     return norm2(dx);
 310 }
 311
 312 /*! \brief Append t_idef structures 1 to nsrc in src to *dest */
 313 static void combine_idef(InteractionDefinitions* dest, gmx::ArrayRef<const thread_work_t> src)
 314 {
 315     for (int ftype = 0; ftype < F_NRE; ftype++)
 316     {
 317         int n = 0;
 318         for (gmx::index s = 1; s < src.ssize(); s++)
 319         {
 320             n += src[s].idef.il[ftype].size();
 321         }
 322         if (n > 0)
 323         {
 324             for (gmx::index s = 1; s < src.ssize(); s++)
 325             {
 326                 dest->il[ftype].append(src[s].idef.il[ftype]);
 327             }
 328
 329             /* Position restraints need an additional treatment */
 330             if (ftype == F_POSRES || ftype == F_FBPOSRES)
 331             {
 332                 int                     nposres = dest->il[ftype].size() / 2;
 333                 std::vector<t_iparams>& iparams_dest =
 334                         (ftype == F_POSRES ? dest->iparams_posres : dest->iparams_fbposres);
 335
 336                 /* Set nposres to the number of original position restraints in dest */
 337                 for (gmx::index s = 1; s < src.ssize(); s++)
 338                 {
 339                     nposres -= src[s].idef.il[ftype].size() / 2;
 340                 }
 341
 342                 for (gmx::index s = 1; s < src.ssize(); s++)
 343                 {
 344                     const std::vector<t_iparams>& iparams_src =
 345                             (ftype == F_POSRES ? src[s].idef.iparams_posres : src[s].idef.iparams_fbposres);
 346                     iparams_dest.insert(iparams_dest.end(), iparams_src.begin(), iparams_src.end());
 347
 348                     /* Correct the indices into iparams_posres */
 349                     for (int i = 0; i < src[s].idef.il[ftype].size() / 2; i++)
 350                     {
 351                         /* Correct the index into iparams_posres */
 352                         dest->il[ftype].iatoms[nposres * 2] = nposres;
 353                         nposres++;
 354                     }
 355                 }
 356                 GMX_RELEASE_ASSERT(
 357                         int(iparams_dest.size()) == nposres,
 358                         "The number of parameters should match the number of restraints");
 359             }
 360         }
 361     }
 362 }
 363
 364 //! Options for assigning interactions for atoms
 365 enum class InteractionConnectivity
 366 {
 367     Intramolecular, //!< Only intra-molecular interactions
 368     Intermolecular  //!< Only inter-molecular interactions
 369 };
 370
 371 /*! \brief Determine whether the local domain has responsibility for
 372  * any of the bonded interactions for local atom \p atomIndex
 373  * and assign those to the local domain.
 374  *
 375  * \returns The total number of bonded interactions for this atom for
 376  * which this domain is responsible.
 377  */
 378 template<InteractionConnectivity interactionConnectivity>
 379 static inline int assignInteractionsForAtom(const int                 atomIndex,
 380                                             const int                 globalAtomIndex,
 381                                             const int                 atomIndexInMolecule,
 382                                             gmx::ArrayRef<const int>  index,
 383                                             gmx::ArrayRef<const int>  rtil,
 384                                             const int                 ind_start,
 385                                             const int                 ind_end,
 386                                             const gmx_ga2la_t&        ga2la,
 387                                             const gmx_domdec_zones_t* zones,
 388                                             const bool                checkDistanceMultiBody,
 389                                             const ivec                rcheck,
 390                                             const bool                checkDistanceTwoBody,
 391                                             const real                cutoffSquared,
 392                                             const t_pbc*              pbc_null,
 393                                             ArrayRef<const RVec>      coordinates,
 394                                             InteractionDefinitions*   idef,
 395                                             const int                 iz,
 396                                             const DDBondedChecking    ddBondedChecking)
 397 {
 398     gmx::ArrayRef<const DDPairInteractionRanges> iZones = zones->iZones;
 399
 400     int numBondedInteractions = 0;
 401
 402     int j = ind_start;
 403     while (j < ind_end)
 404     {
 405         t_iatom tiatoms[1 + MAXATOMLIST];
 406
 407         const int ftype  = rtil[j++];
 408         auto      iatoms = gmx::constArrayRefFromArray(rtil.data() + j, rtil.size() - j);
 409         const int nral   = NRAL(ftype);
 410         if (interaction_function[ftype].flags & IF_VSITE)
 411         {
 412             GMX_ASSERT(interactionConnectivity == InteractionConnectivity::Intramolecular,
 413                        "All vsites should be intramolecular");
 414
 415             /* The vsite construction goes where the vsite itself is */
 416             if (iz == 0)
 417             {
 418                 add_vsite(ga2la,
 419                           index,
 420                           rtil,
 421                           ftype,
 422                           nral,
 423                           TRUE,
 424                           atomIndex,
 425                           globalAtomIndex,
 426                           atomIndexInMolecule,
 427                           iatoms.data(),
 428                           idef);
 429             }
 430         }
 431         else
 432         {
 433             bool bUse = false;
 434
 435             /* Copy the type */
 436             tiatoms[0] = iatoms[0];
 437
 438             if (nral == 1)
 439             {
 440                 GMX_ASSERT(interactionConnectivity == InteractionConnectivity::Intramolecular,
 441                            "All interactions that involve a single atom are intramolecular");
 442
 443                 /* Assign single-body interactions to the home zone.
 444                  * Position restraints are not handled here, but separately.
 445                  */
 446                 if (iz == 0 && !(ftype == F_POSRES || ftype == F_FBPOSRES))
 447                 {
 448                     bUse       = true;
 449                     tiatoms[1] = atomIndex;
 450                 }
 451             }
 452             else if (nral == 2)
 453             {
 454                 /* This is a two-body interaction, we can assign
 455                  * analogous to the non-bonded assignments.
 456                  */
 457                 const int k_gl = (interactionConnectivity == InteractionConnectivity::Intramolecular)
 458                                          ?
 459                                          /* Get the global index using the offset in the molecule */
 460                                          (globalAtomIndex + iatoms[2] - atomIndexInMolecule)
 461                                          : iatoms[2];
 462                 if (const auto* entry = ga2la.find(k_gl))
 463                 {
 464                     int kz = entry->cell;
 465                     if (kz >= zones->n)
 466                     {
 467                         kz -= zones->n;
 468                     }
 469                     /* Check zone interaction assignments */
 470                     bUse = ((iz < iZones.ssize() && iz <= kz && iZones[iz].jZoneRange.isInRange(kz))
 471                             || (kz < iZones.ssize() && iz > kz && iZones[kz].jZoneRange.isInRange(iz)));
 472                     if (bUse)
 473                     {
 474                         GMX_ASSERT(ftype != F_CONSTR || (iz == 0 && kz == 0),
 475                                    "Constraint assigned here should only involve home atoms");
 476
 477                         tiatoms[1] = atomIndex;
 478                         tiatoms[2] = entry->la;
 479                         /* If necessary check the cgcm distance */
 480                         if (checkDistanceTwoBody
 481                             && dd_dist2(pbc_null, coordinates, tiatoms[1], tiatoms[2]) >= cutoffSquared)
 482                         {
 483                             bUse = false;
 484                         }
 485                     }
 486                 }
 487                 else
 488                 {
 489                     bUse = false;
 490                 }
 491             }
 492             else
 493             {
 494                 /* Assign this multi-body bonded interaction to
 495                  * the local domain if we have all the atoms involved
 496                  * (local or communicated) and the minimum zone shift
 497                  * in each dimension is zero, for dimensions
 498                  * with 2 DD cells an extra check may be necessary.
 499                  */
 500                 ivec k_zero, k_plus;
 501                 bUse = true;
 502                 clear_ivec(k_zero);
 503                 clear_ivec(k_plus);
 504                 for (int k = 1; k <= nral && bUse; k++)
 505                 {
 506                     const int k_gl =
 507                             (interactionConnectivity == InteractionConnectivity::Intramolecular)
 508                                     ?
 509                                     /* Get the global index using the offset in the molecule */
 510                                     (globalAtomIndex + iatoms[k] - atomIndexInMolecule)
 511                                     : iatoms[k];
 512                     const auto* entry = ga2la.find(k_gl);
 513                     if (entry == nullptr || entry->cell >= zones->n)
 514                     {
 515                         /* We do not have this atom of this interaction
 516                          * locally, or it comes from more than one cell
 517                          * away.
 518                          */
 519                         bUse = FALSE;
 520                     }
 521                     else
 522                     {
 523                         tiatoms[k] = entry->la;
 524                         for (int d = 0; d < DIM; d++)
 525                         {
 526                             if (zones->shift[entry->cell][d] == 0)
 527                             {
 528                                 k_zero[d] = k;
 529                             }
 530                             else
 531                             {
 532                                 k_plus[d] = k;
 533                             }
 534                         }
 535                     }
 536                 }
 537                 bUse = (bUse && (k_zero[XX] != 0) && (k_zero[YY] != 0) && (k_zero[ZZ] != 0));
 538                 if (checkDistanceMultiBody)
 539                 {
 540                     for (int d = 0; (d < DIM && bUse); d++)
 541                     {
 542                         /* Check if the distance falls within
 543                          * the cut-off to avoid possible multiple
 544                          * assignments of bonded interactions.
 545                          */
 546                         if (rcheck[d] && k_plus[d]
 547                             && dd_dist2(pbc_null, coordinates, tiatoms[k_zero[d]], tiatoms[k_plus[d]])
 548                                        >= cutoffSquared)
 549                         {
 550                             bUse = FALSE;
 551                         }
 552                     }
 553                 }
 554             }
 555             if (bUse)
 556             {
 557                 /* Add this interaction to the local topology */
 558                 idef->il[ftype].push_back(tiatoms[0], nral, tiatoms + 1);
 559                 /* Sum so we can check in global_stat
 560                  * if we have everything.
 561                  */
 562                 if (ddBondedChecking == DDBondedChecking::All
 563                     || !(interaction_function[ftype].flags & IF_LIMZERO))
 564                 {
 565                     numBondedInteractions++;
 566                 }
 567             }
 568         }
 569         j += 1 + nral_rt(ftype);
 570     }
 571
 572     return numBondedInteractions;
 573 }
 574
 575 /*! \brief Determine whether the local domain has responsibility for
 576  * any of the position restraints for local atom \p atomIndex
 577  * and assign those to the local domain.
 578  *
 579  * \returns The total number of bonded interactions for this atom for
 580  * which this domain is responsible.
 581  */
 582 static inline int assignPositionRestraintsForAtom(const int                atomIndex,
 583                                                   const int                mol,
 584                                                   const int                atomIndexInMolecule,
 585                                                   const int                numAtomsInMolecule,
 586                                                   gmx::ArrayRef<const int> rtil,
 587                                                   const int                ind_start,
 588                                                   const int                ind_end,
 589                                                   const gmx_molblock_t&    molb,
 590                                                   const t_iparams*         ip_in,
 591                                                   InteractionDefinitions*  idef)
 592 {
 593     constexpr int nral = 1;
 594
 595     int numBondedInteractions = 0;
 596
 597     int j = ind_start;
 598     while (j < ind_end)
 599     {
 600         const int ftype  = rtil[j++];
 601         auto      iatoms = gmx::constArrayRefFromArray(rtil.data() + j, rtil.size() - j);
 602
 603         if (ftype == F_POSRES || ftype == F_FBPOSRES)
 604         {
 605             std::array<int, 1 + nral> tiatoms = { iatoms[0], atomIndex };
 606             if (ftype == F_POSRES)
 607             {
 608                 add_posres(mol, atomIndexInMolecule, numAtomsInMolecule, molb, tiatoms, ip_in, idef);
 609             }
 610             else
 611             {
 612                 add_fbposres(mol, atomIndexInMolecule, numAtomsInMolecule, molb, tiatoms, ip_in, idef);
 613             }
 614             idef->il[ftype].push_back(tiatoms[0], nral, tiatoms.data() + 1);
 615             numBondedInteractions++;
 616         }
 617         j += 1 + nral_rt(ftype);
 618     }
 619
 620     return numBondedInteractions;
 621 }
 622
 623 /*! \brief This function looks up and assigns bonded interactions for zone iz.
 624  *
 625  * With thread parallelizing each thread acts on a different atom range:
 626  * at_start to at_end.
 627  */
 628 static int make_bondeds_zone(gmx_reverse_top_t*                 rt,
 629                              ArrayRef<const int>                globalAtomIndices,
 630                              const gmx_ga2la_t&                 ga2la,
 631                              const gmx_domdec_zones_t*          zones,
 632                              const std::vector<gmx_molblock_t>& molb,
 633                              const bool                         checkDistanceMultiBody,
 634                              const ivec                         rcheck,
 635                              const bool                         checkDistanceTwoBody,
 636                              const real                         cutoffSquared,
 637                              const t_pbc*                       pbc_null,
 638                              ArrayRef<const RVec>               coordinates,
 639                              const t_iparams*                   ip_in,
 640                              InteractionDefinitions*            idef,
 641                              int                                izone,
 642                              const gmx::Range<int>&             atomRange)
 643 {
 644     int mb                  = 0;
 645     int mt                  = 0;
 646     int mol                 = 0;
 647     int atomIndexInMolecule = 0;
 648
 649     const auto ddBondedChecking = rt->options().ddBondedChecking;
 650
 651     int numBondedInteractions = 0;
 652
 653     for (int i : atomRange)
 654     {
 655         /* Get the global atom number */
 656         const int globalAtomIndex = globalAtomIndices[i];
 657         global_atomnr_to_moltype_ind(
 658                 rt->molblockIndices(), globalAtomIndex, &mb, &mt, &mol, &atomIndexInMolecule);
 659         /* Check all intramolecular interactions assigned to this atom */
 660         gmx::ArrayRef<const int>     index = rt->interactionListForMoleculeType(mt).index;
 661         gmx::ArrayRef<const t_iatom> rtil  = rt->interactionListForMoleculeType(mt).il;
 662
 663         numBondedInteractions += assignInteractionsForAtom<InteractionConnectivity::Intramolecular>(
 664                 i,
 665                 globalAtomIndex,
 666                 atomIndexInMolecule,
 667                 index,
 668                 rtil,
 669                 index[atomIndexInMolecule],
 670                 index[atomIndexInMolecule + 1],
 671                 ga2la,
 672                 zones,
 673                 checkDistanceMultiBody,
 674                 rcheck,
 675                 checkDistanceTwoBody,
 676                 cutoffSquared,
 677                 pbc_null,
 678                 coordinates,
 679                 idef,
 680                 izone,
 681                 ddBondedChecking);
 682
 683         // Assign position restraints, when present, for the home zone
 684         if (izone == 0 && rt->hasPositionRestraints())
 685         {
 686             numBondedInteractions += assignPositionRestraintsForAtom(
 687                     i,
 688                     mol,
 689                     atomIndexInMolecule,
 690                     rt->interactionListForMoleculeType(mt).numAtomsInMolecule,
 691                     rtil,
 692                     index[atomIndexInMolecule],
 693                     index[atomIndexInMolecule + 1],
 694                     molb[mb],
 695                     ip_in,
 696                     idef);
 697         }
 698
 699         if (rt->hasIntermolecularInteractions())
 700         {
 701             /* Check all intermolecular interactions assigned to this atom */
 702             index = rt->interactionListForIntermolecularInteractions().index;
 703             rtil  = rt->interactionListForIntermolecularInteractions().il;
 704
 705             numBondedInteractions += assignInteractionsForAtom<InteractionConnectivity::Intermolecular>(
 706                     i,
 707                     -1, // not used
 708                     -1, // not used
 709                     index,
 710                     rtil,
 711                     index[globalAtomIndex],
 712                     index[globalAtomIndex + 1],
 713                     ga2la,
 714                     zones,
 715                     checkDistanceMultiBody,
 716                     rcheck,
 717                     checkDistanceTwoBody,
 718                     cutoffSquared,
 719                     pbc_null,
 720                     coordinates,
 721                     idef,
 722                     izone,
 723                     ddBondedChecking);
 724         }
 725     }
 726
 727     return numBondedInteractions;
 728 }
 729
 730 /*! \brief Set the exclusion data for i-zone \p iz */
 731 static void make_exclusions_zone(ArrayRef<const int>               globalAtomIndices,
 732                                  const gmx_ga2la_t&                ga2la,
 733                                  gmx_domdec_zones_t*               zones,
 734                                  ArrayRef<const MolblockIndices>   molblockIndices,
 735                                  const std::vector<gmx_moltype_t>& moltype,
 736                                  const int*                        cginfo,
 737                                  ListOfLists<int>*                 lexcls,
 738                                  int                               iz,
 739                                  int                               at_start,
 740                                  int                               at_end,
 741                                  const gmx::ArrayRef<const int>    intermolecularExclusionGroup)
 742 {
 743     const auto& jAtomRange = zones->iZones[iz].jAtomRange;
 744
 745     const gmx::index oldNumLists = lexcls->ssize();
 746
 747     std::vector<int> exclusionsForAtom;
 748     for (int at = at_start; at < at_end; at++)
 749     {
 750         exclusionsForAtom.clear();
 751
 752         if (GET_CGINFO_EXCL_INTER(cginfo[at]))
 753         {
 754             int mb    = 0;
 755             int mt    = 0;
 756             int mol   = 0;
 757             int a_mol = 0;
 758
 759             /* Copy the exclusions from the global top */
 760             int a_gl = globalAtomIndices[at];
 761             global_atomnr_to_moltype_ind(molblockIndices, a_gl, &mb, &mt, &mol, &a_mol);
 762             const auto excls = moltype[mt].excls[a_mol];
 763             for (const int aj_mol : excls)
 764             {
 765                 if (const auto* jEntry = ga2la.find(a_gl + aj_mol - a_mol))
 766                 {
 767                     /* This check is not necessary, but it can reduce
 768                      * the number of exclusions in the list, which in turn
 769                      * can speed up the pair list construction a bit.
 770                      */
 771                     if (jAtomRange.isInRange(jEntry->la))
 772                     {
 773                         exclusionsForAtom.push_back(jEntry->la);
 774                     }
 775                 }
 776             }
 777         }
 778
 779         bool isExcludedAtom = !intermolecularExclusionGroup.empty()
 780                               && std::find(intermolecularExclusionGroup.begin(),
 781                                            intermolecularExclusionGroup.end(),
 782                                            globalAtomIndices[at])
 783                                          != intermolecularExclusionGroup.end();
 784
 785         if (isExcludedAtom)
 786         {
 787             for (int qmAtomGlobalIndex : intermolecularExclusionGroup)
 788             {
 789                 if (const auto* entry = ga2la.find(qmAtomGlobalIndex))
 790                 {
 791                     exclusionsForAtom.push_back(entry->la);
 792                 }
 793             }
 794         }
 795
 796         /* Append the exclusions for this atom to the topology */
 797         lexcls->pushBack(exclusionsForAtom);
 798     }
 799
 800     GMX_RELEASE_ASSERT(
 801             lexcls->ssize() - oldNumLists == at_end - at_start,
 802             "The number of exclusion list should match the number of atoms in the range");
 803 }
 804
 805 /*! \brief Generate and store all required local bonded interactions in \p idef and local exclusions in \p lexcls
 806  *
 807  * \returns Total count of bonded interactions in the local topology on this domain */
 808 static int make_local_bondeds_excls(gmx_domdec_t*           dd,
 809                                     gmx_domdec_zones_t*     zones,
 810                                     const gmx_mtop_t&       mtop,
 811                                     const int*              cginfo,
 812                                     const bool              checkDistanceMultiBody,
 813                                     const ivec              rcheck,
 814                                     const gmx_bool          checkDistanceTwoBody,
 815                                     const real              cutoff,
 816                                     const t_pbc*            pbc_null,
 817                                     ArrayRef<const RVec>    coordinates,
 818                                     InteractionDefinitions* idef,
 819                                     ListOfLists<int>*       lexcls)
 820 {
 821     int nzone_bondeds = 0;
 822
 823     if (dd->reverse_top->hasInterAtomicInteractions())
 824     {
 825         nzone_bondeds = zones->n;
 826     }
 827     else
 828     {
 829         /* Only single charge group (or atom) molecules, so interactions don't
 830          * cross zone boundaries and we only need to assign in the home zone.
 831          */
 832         nzone_bondeds = 1;
 833     }
 834
 835     /* We only use exclusions from i-zones to i- and j-zones */
 836     const int numIZonesForExclusions = (dd->haveExclusions ? zones->iZones.size() : 0);
 837
 838     gmx_reverse_top_t* rt = dd->reverse_top.get();
 839
 840     const real cutoffSquared = gmx::square(cutoff);
 841
 842     /* Clear the counts */
 843     idef->clear();
 844     int numBondedInteractions = 0;
 845
 846     lexcls->clear();
 847
 848     for (int izone = 0; izone < nzone_bondeds; izone++)
 849     {
 850         const int cg0 = zones->cg_range[izone];
 851         const int cg1 = zones->cg_range[izone + 1];
 852
 853         gmx::ArrayRef<thread_work_t> threadWorkObjects = rt->threadWorkObjects();
 854         const int                    numThreads        = threadWorkObjects.size();
 855 #pragma omp parallel for num_threads(numThreads) schedule(static)
 856         for (int thread = 0; thread < numThreads; thread++)
 857         {
 858             try
 859             {
 860                 InteractionDefinitions* idef_t = nullptr;
 861
 862                 int cg0t = cg0 + ((cg1 - cg0) * thread) / numThreads;
 863                 int cg1t = cg0 + ((cg1 - cg0) * (thread + 1)) / numThreads;
 864
 865                 if (thread == 0)
 866                 {
 867                     idef_t = idef;
 868                 }
 869                 else
 870                 {
 871                     idef_t = &threadWorkObjects[thread].idef;
 872                     idef_t->clear();
 873                 }
 874
 875                 threadWorkObjects[thread].numBondedInteractions =
 876                         make_bondeds_zone(rt,
 877                                           dd->globalAtomIndices,
 878                                           *dd->ga2la,
 879                                           zones,
 880                                           mtop.molblock,
 881                                           checkDistanceMultiBody,
 882                                           rcheck,
 883                                           checkDistanceTwoBody,
 884                                           cutoffSquared,
 885                                           pbc_null,
 886                                           coordinates,
 887                                           idef->iparams.data(),
 888                                           idef_t,
 889                                           izone,
 890                                           gmx::Range<int>(cg0t, cg1t));
 891
 892                 if (izone < numIZonesForExclusions)
 893                 {
 894                     ListOfLists<int>* excl_t = nullptr;
 895                     if (thread == 0)
 896                     {
 897                         // Thread 0 stores exclusions directly in the final storage
 898                         excl_t = lexcls;
 899                     }
 900                     else
 901                     {
 902                         // Threads > 0 store in temporary storage, starting at list index 0
 903                         excl_t = &threadWorkObjects[thread].excl;
 904                         excl_t->clear();
 905                     }
 906
 907                     /* No charge groups and no distance check required */
 908                     make_exclusions_zone(dd->globalAtomIndices,
 909                                          *dd->ga2la,
 910                                          zones,
 911                                          rt->molblockIndices(),
 912                                          mtop.moltype,
 913                                          cginfo,
 914                                          excl_t,
 915                                          izone,
 916                                          cg0t,
 917                                          cg1t,
 918                                          mtop.intermolecularExclusionGroup);
 919                 }
 920             }
 921             GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
 922         }
 923
 924         if (threadWorkObjects.size() > 1)
 925         {
 926             combine_idef(idef, threadWorkObjects);
 927         }
 928
 929         for (const thread_work_t& th_work : threadWorkObjects)
 930         {
 931             numBondedInteractions += th_work.numBondedInteractions;
 932         }
 933
 934         if (izone < numIZonesForExclusions)
 935         {
 936             for (std::size_t th = 1; th < threadWorkObjects.size(); th++)
 937             {
 938                 lexcls->appendListOfLists(threadWorkObjects[th].excl);
 939             }
 940         }
 941     }
 942
 943     if (debug)
 944     {
 945         fprintf(debug, "We have %d exclusions\n", lexcls->numElements());
 946     }
 947
 948     return numBondedInteractions;
 949 }
 950
 951 int dd_make_local_top(gmx_domdec_t*        dd,
 952                       gmx_domdec_zones_t*  zones,
 953                       int                  npbcdim,
 954                       matrix               box,
 955                       rvec                 cellsize_min,
 956                       const ivec           npulse,
 957                       t_forcerec*          fr,
 958                       ArrayRef<const RVec> coordinates,
 959                       const gmx_mtop_t&    mtop,
 960                       gmx_localtop_t*      ltop)
 961 {
 962     real  rc = -1;
 963     ivec  rcheck;
 964     t_pbc pbc, *pbc_null = nullptr;
 965
 966     if (debug)
 967     {
 968         fprintf(debug, "Making local topology\n");
 969     }
 970
 971     bool checkDistanceMultiBody = false;
 972     bool checkDistanceTwoBody   = false;
 973
 974     if (dd->reverse_top->hasInterAtomicInteractions())
 975     {
 976         /* We need to check to which cell bondeds should be assigned */
 977         rc = dd_cutoff_twobody(dd);
 978         if (debug)
 979         {
 980             fprintf(debug, "Two-body bonded cut-off distance is %g\n", rc);
 981         }
 982
 983         /* Should we check distances when assigning bonded interactions? */
 984         for (int d = 0; d < DIM; d++)
 985         {
 986             rcheck[d] = FALSE;
 987             /* Only need to check for dimensions where the part of the box
 988              * that is not communicated is smaller than the cut-off.
 989              */
 990             if (d < npbcdim && dd->numCells[d] > 1
 991                 && (dd->numCells[d] - npulse[d]) * cellsize_min[d] < 2 * rc)
 992             {
 993                 if (dd->numCells[d] == 2)
 994                 {
 995                     rcheck[d]              = TRUE;
 996                     checkDistanceMultiBody = TRUE;
 997                 }
 998                 /* Check for interactions between two atoms,
 999                  * where we can allow interactions up to the cut-off,
1000                  * instead of up to the smallest cell dimension.
1001                  */
1002                 checkDistanceTwoBody = TRUE;
1003             }
1004             if (debug)
1005             {
1006                 fprintf(debug,
1007                         "dim %d cellmin %f bonded rcheck[%d] = %d, checkDistanceTwoBody = %s\n",
1008                         d,
1009                         cellsize_min[d],
1010                         d,
1011                         rcheck[d],
1012                         gmx::boolToString(checkDistanceTwoBody));
1013             }
1014         }
1015         if (checkDistanceMultiBody || checkDistanceTwoBody)
1016         {
1017             if (fr->bMolPBC)
1018             {
1019                 pbc_null = set_pbc_dd(&pbc, fr->pbcType, dd->numCells, TRUE, box);
1020             }
1021             else
1022             {
1023                 pbc_null = nullptr;
1024             }
1025         }
1026     }
1027
1028     int numBondedInteractionsToReduce = make_local_bondeds_excls(dd,
1029                                                                  zones,
1030                                                                  mtop,
1031                                                                  fr->cginfo.data(),
1032                                                                  checkDistanceMultiBody,
1033                                                                  rcheck,
1034                                                                  checkDistanceTwoBody,
1035                                                                  rc,
1036                                                                  pbc_null,
1037                                                                  coordinates,
1038                                                                  &ltop->idef,
1039                                                                  &ltop->excls);
1040
1041     /* The ilist is not sorted yet,
1042      * we can only do this when we have the charge arrays.
1043      */
1044     ltop->idef.ilsort = ilsortUNKNOWN;
1045
1046     return numBondedInteractionsToReduce;
1047 }
1048
1049 void dd_sort_local_top(const gmx_domdec_t& dd, const t_mdatoms* mdatoms, gmx_localtop_t* ltop)
1050 {
1051     if (dd.reverse_top->doSorting())
1052     {
1053         gmx_sort_ilist_fe(&ltop->idef, mdatoms->chargeA, mdatoms->chargeB);
1054     }
1055     else
1056     {
1057         ltop->idef.ilsort = ilsortNO_FE;
1058     }
1059 }
1060
1061 void dd_init_local_state(const gmx_domdec_t& dd, const t_state* state_global, t_state* state_local)
1062 {
1063     int buf[NITEM_DD_INIT_LOCAL_STATE];
1064
1065     if (DDMASTER(dd))
1066     {
1067         buf[0] = state_global->flags;
1068         buf[1] = state_global->ngtc;
1069         buf[2] = state_global->nnhpres;
1070         buf[3] = state_global->nhchainlength;
1071         buf[4] = state_global->dfhist ? state_global->dfhist->nlambda : 0;
1072     }
1073     dd_bcast(&dd, NITEM_DD_INIT_LOCAL_STATE * sizeof(int), buf);
1074
1075     init_gtc_state(state_local, buf[1], buf[2], buf[3]);
1076     init_dfhist_state(state_local, buf[4]);
1077     state_local->flags = buf[0];
1078 }
1079
1080 /*! \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 */
1081 static void check_link(t_blocka* link, int cg_gl, int cg_gl_j)
1082 {
1083     bool bFound = false;
1084     for (int k = link->index[cg_gl]; k < link->index[cg_gl + 1]; k++)
1085     {
1086         GMX_RELEASE_ASSERT(link->a, "Inconsistent NULL pointer while making charge-group links");
1087         if (link->a[k] == cg_gl_j)
1088         {
1089             bFound = TRUE;
1090         }
1091     }
1092     if (!bFound)
1093     {
1094         GMX_RELEASE_ASSERT(link->a || link->index[cg_gl + 1] + 1 > link->nalloc_a,
1095                            "Inconsistent allocation of link");
1096         /* Add this charge group link */
1097         if (link->index[cg_gl + 1] + 1 > link->nalloc_a)
1098         {
1099             link->nalloc_a = over_alloc_large(link->index[cg_gl + 1] + 1);
1100             srenew(link->a, link->nalloc_a);
1101         }
1102         link->a[link->index[cg_gl + 1]] = cg_gl_j;
1103         link->index[cg_gl + 1]++;
1104     }
1105 }
1106
1107 void makeBondedLinks(gmx_domdec_t* dd, const gmx_mtop_t& mtop, gmx::ArrayRef<cginfo_mb_t> cginfo_mb)
1108 {
1109
1110     if (!dd->comm->systemInfo.filterBondedCommunication)
1111     {
1112         /* Only communicate atoms based on cut-off */
1113         dd->comm->bondedLinks = nullptr;
1114         return;
1115     }
1116
1117     t_blocka* link = nullptr;
1118
1119     /* For each atom make a list of other atoms in the system
1120      * that a linked to it via bonded interactions
1121      * which are also stored in reverse_top.
1122      */
1123
1124     reverse_ilist_t ril_intermol;
1125     if (mtop.bIntermolecularInteractions)
1126     {
1127         t_atoms atoms;
1128
1129         atoms.nr   = mtop.natoms;
1130         atoms.atom = nullptr;
1131
1132         GMX_RELEASE_ASSERT(mtop.intermolecular_ilist,
1133                            "We should have an ilist when intermolecular interactions are on");
1134
1135         ReverseTopOptions rtOptions(DDBondedChecking::ExcludeZeroLimit);
1136         make_reverse_ilist(
1137                 *mtop.intermolecular_ilist, &atoms, rtOptions, AtomLinkRule::AllAtomsInBondeds, &ril_intermol);
1138     }
1139
1140     snew(link, 1);
1141     snew(link->index, mtop.natoms + 1);
1142     link->nalloc_a = 0;
1143     link->a        = nullptr;
1144
1145     link->index[0] = 0;
1146     int cg_offset  = 0;
1147     int ncgi       = 0;
1148     for (size_t mb = 0; mb < mtop.molblock.size(); mb++)
1149     {
1150         const gmx_molblock_t& molb = mtop.molblock[mb];
1151         if (molb.nmol == 0)
1152         {
1153             continue;
1154         }
1155         const gmx_moltype_t& molt = mtop.moltype[molb.type];
1156         /* Make a reverse ilist in which the interactions are linked
1157          * to all atoms, not only the first atom as in gmx_reverse_top.
1158          * The constraints are discarded here.
1159          */
1160         ReverseTopOptions rtOptions(DDBondedChecking::ExcludeZeroLimit);
1161         reverse_ilist_t   ril;
1162         make_reverse_ilist(molt.ilist, &molt.atoms, rtOptions, AtomLinkRule::AllAtomsInBondeds, &ril);
1163
1164         cginfo_mb_t* cgi_mb = &cginfo_mb[mb];
1165
1166         int mol = 0;
1167         for (mol = 0; mol < (mtop.bIntermolecularInteractions ? molb.nmol : 1); mol++)
1168         {
1169             for (int a = 0; a < molt.atoms.nr; a++)
1170             {
1171                 int cg_gl              = cg_offset + a;
1172                 link->index[cg_gl + 1] = link->index[cg_gl];
1173                 int i                  = ril.index[a];
1174                 while (i < ril.index[a + 1])
1175                 {
1176                     int ftype = ril.il[i++];
1177                     int nral  = NRAL(ftype);
1178                     /* Skip the ifunc index */
1179                     i++;
1180                     for (int j = 0; j < nral; j++)
1181                     {
1182                         int aj = ril.il[i + j];
1183                         if (aj != a)
1184                         {
1185                             check_link(link, cg_gl, cg_offset + aj);
1186                         }
1187                     }
1188                     i += nral_rt(ftype);
1189                 }
1190
1191                 if (mtop.bIntermolecularInteractions)
1192                 {
1193                     int i = ril_intermol.index[cg_gl];
1194                     while (i < ril_intermol.index[cg_gl + 1])
1195                     {
1196                         int ftype = ril_intermol.il[i++];
1197                         int nral  = NRAL(ftype);
1198                         /* Skip the ifunc index */
1199                         i++;
1200                         for (int j = 0; j < nral; j++)
1201                         {
1202                             /* Here we assume we have no charge groups;
1203                              * this has been checked above.
1204                              */
1205                             int aj = ril_intermol.il[i + j];
1206                             check_link(link, cg_gl, aj);
1207                         }
1208                         i += nral_rt(ftype);
1209                     }
1210                 }
1211                 if (link->index[cg_gl + 1] - link->index[cg_gl] > 0)
1212                 {
1213                     SET_CGINFO_BOND_INTER(cgi_mb->cginfo[a]);
1214                     ncgi++;
1215                 }
1216             }
1217
1218             cg_offset += molt.atoms.nr;
1219         }
1220         int nlink_mol = link->index[cg_offset] - link->index[cg_offset - molt.atoms.nr];
1221
1222         if (debug)
1223         {
1224             fprintf(debug,
1225                     "molecule type '%s' %d atoms has %d atom links through bonded interac.\n",
1226                     *molt.name,
1227                     molt.atoms.nr,
1228                     nlink_mol);
1229         }
1230
1231         if (molb.nmol > mol)
1232         {
1233             /* Copy the data for the rest of the molecules in this block */
1234             link->nalloc_a += (molb.nmol - mol) * nlink_mol;
1235             srenew(link->a, link->nalloc_a);
1236             for (; mol < molb.nmol; mol++)
1237             {
1238                 for (int a = 0; a < molt.atoms.nr; a++)
1239                 {
1240                     int cg_gl              = cg_offset + a;
1241                     link->index[cg_gl + 1] = link->index[cg_gl + 1 - molt.atoms.nr] + nlink_mol;
1242                     for (int j = link->index[cg_gl]; j < link->index[cg_gl + 1]; j++)
1243                     {
1244                         link->a[j] = link->a[j - nlink_mol] + molt.atoms.nr;
1245                     }
1246                     if (link->index[cg_gl + 1] - link->index[cg_gl] > 0
1247                         && cg_gl - cgi_mb->cg_start < cgi_mb->cg_mod)
1248                     {
1249                         SET_CGINFO_BOND_INTER(cgi_mb->cginfo[cg_gl - cgi_mb->cg_start]);
1250                         ncgi++;
1251                     }
1252                 }
1253                 cg_offset += molt.atoms.nr;
1254             }
1255         }
1256     }
1257
1258     if (debug)
1259     {
1260         fprintf(debug, "Of the %d atoms %d are linked via bonded interactions\n", mtop.natoms, ncgi);
1261     }
1262
1263     dd->comm->bondedLinks = link;
1264 }
1265
1266 typedef struct
1267 {
1268     real r2;
1269     int  ftype;
1270     int  a1;
1271     int  a2;
1272 } bonded_distance_t;
1273
1274 /*! \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 */
1275 static void update_max_bonded_distance(real r2, int ftype, int a1, int a2, bonded_distance_t* bd)
1276 {
1277     if (r2 > bd->r2)
1278     {
1279         bd->r2    = r2;
1280         bd->ftype = ftype;
1281         bd->a1    = a1;
1282         bd->a2    = a2;
1283     }
1284 }
1285
1286 /*! \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 */
1287 static void bonded_cg_distance_mol(const gmx_moltype_t*   molt,
1288                                    const DDBondedChecking ddBondedChecking,
1289                                    gmx_bool               bExcl,
1290                                    ArrayRef<const RVec>   x,
1291                                    bonded_distance_t*     bd_2b,
1292                                    bonded_distance_t*     bd_mb)
1293 {
1294     const ReverseTopOptions rtOptions(ddBondedChecking);
1295
1296     for (int ftype = 0; ftype < F_NRE; ftype++)
1297     {
1298         if (dd_check_ftype(ftype, rtOptions))
1299         {
1300             const auto& il   = molt->ilist[ftype];
1301             int         nral = NRAL(ftype);
1302             if (nral > 1)
1303             {
1304                 for (int i = 0; i < il.size(); i += 1 + nral)
1305                 {
1306                     for (int ai = 0; ai < nral; ai++)
1307                     {
1308                         int atomI = il.iatoms[i + 1 + ai];
1309                         for (int aj = ai + 1; aj < nral; aj++)
1310                         {
1311                             int atomJ = il.iatoms[i + 1 + aj];
1312                             if (atomI != atomJ)
1313                             {
1314                                 real rij2 = distance2(x[atomI], x[atomJ]);
1315
1316                                 update_max_bonded_distance(
1317                                         rij2, ftype, atomI, atomJ, (nral == 2) ? bd_2b : bd_mb);
1318                             }
1319                         }
1320                     }
1321                 }
1322             }
1323         }
1324     }
1325     if (bExcl)
1326     {
1327         const auto& excls = molt->excls;
1328         for (gmx::index ai = 0; ai < excls.ssize(); ai++)
1329         {
1330             for (const int aj : excls[ai])
1331             {
1332                 if (ai != aj)
1333                 {
1334                     real rij2 = distance2(x[ai], x[aj]);
1335
1336                     /* There is no function type for exclusions, use -1 */
1337                     update_max_bonded_distance(rij2, -1, ai, aj, bd_2b);
1338                 }
1339             }
1340         }
1341     }
1342 }
1343
1344 /*! \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 */
1345 static void bonded_distance_intermol(const InteractionLists& ilists_intermol,
1346                                      const DDBondedChecking  ddBondedChecking,
1347                                      ArrayRef<const RVec>    x,
1348                                      PbcType                 pbcType,
1349                                      const matrix            box,
1350                                      bonded_distance_t*      bd_2b,
1351                                      bonded_distance_t*      bd_mb)
1352 {
1353     t_pbc pbc;
1354
1355     set_pbc(&pbc, pbcType, box);
1356
1357     const ReverseTopOptions rtOptions(ddBondedChecking);
1358
1359     for (int ftype = 0; ftype < F_NRE; ftype++)
1360     {
1361         if (dd_check_ftype(ftype, rtOptions))
1362         {
1363             const auto& il   = ilists_intermol[ftype];
1364             int         nral = NRAL(ftype);
1365
1366             /* No nral>1 check here, since intermol interactions always
1367              * have nral>=2 (and the code is also correct for nral=1).
1368              */
1369             for (int i = 0; i < il.size(); i += 1 + nral)
1370             {
1371                 for (int ai = 0; ai < nral; ai++)
1372                 {
1373                     int atom_i = il.iatoms[i + 1 + ai];
1374
1375                     for (int aj = ai + 1; aj < nral; aj++)
1376                     {
1377                         rvec dx;
1378
1379                         int atom_j = il.iatoms[i + 1 + aj];
1380
1381                         pbc_dx(&pbc, x[atom_i], x[atom_j], dx);
1382
1383                         const real rij2 = norm2(dx);
1384
1385                         update_max_bonded_distance(rij2, ftype, atom_i, atom_j, (nral == 2) ? bd_2b : bd_mb);
1386                     }
1387                 }
1388             }
1389         }
1390     }
1391 }
1392
1393 //! Returns whether \p molt has at least one virtual site
1394 static bool moltypeHasVsite(const gmx_moltype_t& molt)
1395 {
1396     bool hasVsite = false;
1397     for (int i = 0; i < F_NRE; i++)
1398     {
1399         if ((interaction_function[i].flags & IF_VSITE) && !molt.ilist[i].empty())
1400         {
1401             hasVsite = true;
1402         }
1403     }
1404
1405     return hasVsite;
1406 }
1407
1408 //! Returns coordinates not broken over PBC for a molecule
1409 static void getWholeMoleculeCoordinates(const gmx_moltype_t*  molt,
1410                                         const gmx_ffparams_t* ffparams,
1411                                         PbcType               pbcType,
1412                                         t_graph*              graph,
1413                                         const matrix          box,
1414                                         ArrayRef<const RVec>  x,
1415                                         ArrayRef<RVec>        xs)
1416 {
1417     if (pbcType != PbcType::No)
1418     {
1419         mk_mshift(nullptr, graph, pbcType, box, as_rvec_array(x.data()));
1420
1421         shift_x(graph, box, as_rvec_array(x.data()), as_rvec_array(xs.data()));
1422         /* By doing an extra mk_mshift the molecules that are broken
1423          * because they were e.g. imported from another software
1424          * will be made whole again. Such are the healing powers
1425          * of GROMACS.
1426          */
1427         mk_mshift(nullptr, graph, pbcType, box, as_rvec_array(xs.data()));
1428     }
1429     else
1430     {
1431         /* We copy the coordinates so the original coordinates remain
1432          * unchanged, just to be 100% sure that we do not affect
1433          * binary reproducibility of simulations.
1434          */
1435         for (int i = 0; i < molt->atoms.nr; i++)
1436         {
1437             copy_rvec(x[i], xs[i]);
1438         }
1439     }
1440
1441     if (moltypeHasVsite(*molt))
1442     {
1443         gmx::constructVirtualSites(xs, ffparams->iparams, molt->ilist);
1444     }
1445 }
1446
1447 void dd_bonded_cg_distance(const gmx::MDLogger&   mdlog,
1448                            const gmx_mtop_t&      mtop,
1449                            const t_inputrec&      inputrec,
1450                            ArrayRef<const RVec>   x,
1451                            const matrix           box,
1452                            const DDBondedChecking ddBondedChecking,
1453                            real*                  r_2b,
1454                            real*                  r_mb)
1455 {
1456     bonded_distance_t bd_2b = { 0, -1, -1, -1 };
1457     bonded_distance_t bd_mb = { 0, -1, -1, -1 };
1458
1459     bool bExclRequired = inputrecExclForces(&inputrec);
1460
1461     *r_2b         = 0;
1462     *r_mb         = 0;
1463     int at_offset = 0;
1464     for (const gmx_molblock_t& molb : mtop.molblock)
1465     {
1466         const gmx_moltype_t& molt = mtop.moltype[molb.type];
1467         if (molt.atoms.nr == 1 || molb.nmol == 0)
1468         {
1469             at_offset += molb.nmol * molt.atoms.nr;
1470         }
1471         else
1472         {
1473             t_graph graph;
1474             if (inputrec.pbcType != PbcType::No)
1475             {
1476                 graph = mk_graph_moltype(molt);
1477             }
1478
1479             std::vector<RVec> xs(molt.atoms.nr);
1480             for (int mol = 0; mol < molb.nmol; mol++)
1481             {
1482                 getWholeMoleculeCoordinates(&molt,
1483                                             &mtop.ffparams,
1484                                             inputrec.pbcType,
1485                                             &graph,
1486                                             box,
1487                                             x.subArray(at_offset, molt.atoms.nr),
1488                                             xs);
1489
1490                 bonded_distance_t bd_mol_2b = { 0, -1, -1, -1 };
1491                 bonded_distance_t bd_mol_mb = { 0, -1, -1, -1 };
1492
1493                 bonded_cg_distance_mol(&molt, ddBondedChecking, bExclRequired, xs, &bd_mol_2b, &bd_mol_mb);
1494
1495                 /* Process the mol data adding the atom index offset */
1496                 update_max_bonded_distance(bd_mol_2b.r2,
1497                                            bd_mol_2b.ftype,
1498                                            at_offset + bd_mol_2b.a1,
1499                                            at_offset + bd_mol_2b.a2,
1500                                            &bd_2b);
1501                 update_max_bonded_distance(bd_mol_mb.r2,
1502                                            bd_mol_mb.ftype,
1503                                            at_offset + bd_mol_mb.a1,
1504                                            at_offset + bd_mol_mb.a2,
1505                                            &bd_mb);
1506
1507                 at_offset += molt.atoms.nr;
1508             }
1509         }
1510     }
1511
1512     if (mtop.bIntermolecularInteractions)
1513     {
1514         GMX_RELEASE_ASSERT(mtop.intermolecular_ilist,
1515                            "We should have an ilist when intermolecular interactions are on");
1516
1517         bonded_distance_intermol(
1518                 *mtop.intermolecular_ilist, ddBondedChecking, x, inputrec.pbcType, box, &bd_2b, &bd_mb);
1519     }
1520
1521     *r_2b = sqrt(bd_2b.r2);
1522     *r_mb = sqrt(bd_mb.r2);
1523
1524     if (*r_2b > 0 || *r_mb > 0)
1525     {
1526         GMX_LOG(mdlog.info).appendText("Initial maximum distances in bonded interactions:");
1527         if (*r_2b > 0)
1528         {
1529             GMX_LOG(mdlog.info)
1530                     .appendTextFormatted(
1531                             "    two-body bonded interactions: %5.3f nm, %s, atoms %d %d",
1532                             *r_2b,
1533                             (bd_2b.ftype >= 0) ? interaction_function[bd_2b.ftype].longname : "Exclusion",
1534                             bd_2b.a1 + 1,
1535                             bd_2b.a2 + 1);
1536         }
1537         if (*r_mb > 0)
1538         {
1539             GMX_LOG(mdlog.info)
1540                     .appendTextFormatted(
1541                             "  multi-body bonded interactions: %5.3f nm, %s, atoms %d %d",
1542                             *r_mb,
1543                             interaction_function[bd_mb.ftype].longname,
1544                             bd_mb.a1 + 1,
1545                             bd_mb.a2 + 1);
1546         }
1547     }
1548 }