src/gromacs/domdec/domdec_topology.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2006,2007,2008,2009,2010,2011,2012,2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35
  36 /*! \internal \file
  37  *
  38  * \brief This file defines functions used by the domdec module
  39  * while managing the construction, use and error checking for
  40  * topologies local to a DD rank.
  41  *
  42  * \author Berk Hess <hess@kth.se>
  43  * \ingroup module_domdec
  44  */
  45
  46 #include "gmxpre.h"
  47
  48 #include <assert.h>
  49 #include <stdlib.h>
  50 #include <string.h>
  51
  52 #include <cassert>
  53
  54 #include <algorithm>
  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/chargegroup.h"
  61 #include "gromacs/gmxlib/network.h"
  62 #include "gromacs/math/vec.h"
  63 #include "gromacs/mdlib/force.h"
  64 #include "gromacs/mdlib/forcerec.h"
  65 #include "gromacs/mdlib/gmx_omp_nthreads.h"
  66 #include "gromacs/mdlib/vsite.h"
  67 #include "gromacs/mdtypes/commrec.h"
  68 #include "gromacs/mdtypes/inputrec.h"
  69 #include "gromacs/mdtypes/md_enums.h"
  70 #include "gromacs/mdtypes/mdatom.h"
  71 #include "gromacs/mdtypes/state.h"
  72 #include "gromacs/pbcutil/mshift.h"
  73 #include "gromacs/pbcutil/pbc.h"
  74 #include "gromacs/topology/ifunc.h"
  75 #include "gromacs/topology/mtop_util.h"
  76 #include "gromacs/topology/topsort.h"
  77 #include "gromacs/utility/cstringutil.h"
  78 #include "gromacs/utility/exceptions.h"
  79 #include "gromacs/utility/fatalerror.h"
  80 #include "gromacs/utility/gmxassert.h"
  81 #include "gromacs/utility/smalloc.h"
  82 #include "gromacs/utility/stringutil.h"
  83
  84 #include "domdec_constraints.h"
  85 #include "domdec_internal.h"
  86 #include "domdec_vsite.h"
  87
  88 /*! \brief The number of integer item in the local state, used for broadcasting of the state */
  89 #define NITEM_DD_INIT_LOCAL_STATE 5
  90
  91 typedef struct {
  92     int  *index;  /* Index for each atom into il                  */
  93     int  *il;     /* ftype|type|a0|...|an|ftype|...               */
  94 } reverse_ilist_t;
  95
  96 typedef struct {
  97     int  a_start;
  98     int  a_end;
  99     int  natoms_mol;
 100     int  type;
 101 } molblock_ind_t;
 102
 103 /*! \brief Struct for thread local work data for local topology generation */
 104 typedef struct {
 105     t_idef     idef;             /**< Partial local topology */
 106     int      **vsite_pbc;        /**< vsite PBC structure */
 107     int       *vsite_pbc_nalloc; /**< Allocation sizes for vsite_pbc */
 108     int        nbonded;          /**< The number of bondeds in this struct */
 109     t_blocka   excl;             /**< List of exclusions */
 110     int        excl_count;       /**< The total exclusion count for \p excl */
 111 } thread_work_t;
 112
 113 /*! \brief Struct for the reverse topology: links bonded interactions to atomsx */
 114 struct gmx_reverse_top_t
 115 {
 116     //! @cond Doxygen_Suppress
 117     gmx_bool         bExclRequired;               /**< Do we require all exclusions to be assigned? */
 118     int              n_excl_at_max;               /**< The maximum number of exclusions one atom can have */
 119     gmx_bool         bConstr;                     /**< Are there constraints in this revserse top?  */
 120     gmx_bool         bSettle;                     /**< Are there settles in this revserse top?  */
 121     gmx_bool         bBCheck;                     /**< All bonded interactions have to be assigned? */
 122     gmx_bool         bInterCGInteractions;        /**< Are there bondeds/exclusions between charge-groups? */
 123     reverse_ilist_t *ril_mt;                      /**< Reverse ilist for all moltypes      */
 124     int              ril_mt_tot_size;             /**< The size of ril_mt[?].index summed over all entries */
 125     int              ilsort;                      /**< The sorting state of bondeds for free energy */
 126     molblock_ind_t  *mbi;                         /**< molblock to global atom index for quick lookup of molblocks on atom index */
 127     int              nmolblock;                   /**< The number of molblocks, size of \p mbi */
 128
 129     gmx_bool         bIntermolecularInteractions; /**< Do we have intermolecular interactions? */
 130     reverse_ilist_t  ril_intermol;                /**< Intermolecular reverse ilist */
 131
 132     /* Work data structures for multi-threading */
 133     int            nthread;           /**< The number of threads to be used */
 134     thread_work_t *th_work;           /**< Thread work array for local topology generation */
 135     //! @endcond
 136 };
 137
 138 /*! \brief Returns the number of atom entries for il in gmx_reverse_top_t */
 139 static int nral_rt(int ftype)
 140 {
 141     int nral;
 142
 143     nral = NRAL(ftype);
 144     if (interaction_function[ftype].flags & IF_VSITE)
 145     {
 146         /* With vsites the reverse topology contains
 147          * two extra entries for PBC.
 148          */
 149         nral += 2;
 150     }
 151
 152     return nral;
 153 }
 154
 155 /*! \brief Return whether interactions of type \p ftype need to be assigned exactly once */
 156 static gmx_bool dd_check_ftype(int ftype, gmx_bool bBCheck,
 157                                gmx_bool bConstr, gmx_bool bSettle)
 158 {
 159     return (((interaction_function[ftype].flags & IF_BOND) &&
 160              !(interaction_function[ftype].flags & IF_VSITE) &&
 161              (bBCheck || !(interaction_function[ftype].flags & IF_LIMZERO))) ||
 162             (bConstr && (ftype == F_CONSTR || ftype == F_CONSTRNC)) ||
 163             (bSettle && ftype == F_SETTLE));
 164 }
 165
 166 /*! \brief Print a header on error messages */
 167 static void print_error_header(FILE *fplog, const char *moltypename, int nprint)
 168 {
 169     fprintf(fplog, "\nMolecule type '%s'\n", moltypename);
 170     fprintf(stderr, "\nMolecule type '%s'\n", moltypename);
 171     fprintf(fplog,
 172             "the first %d missing interactions, except for exclusions:\n",
 173             nprint);
 174     fprintf(stderr,
 175             "the first %d missing interactions, except for exclusions:\n",
 176             nprint);
 177 }
 178
 179 /*! \brief Help print error output when interactions are missing */
 180 static void print_missing_interactions_mb(FILE *fplog, t_commrec *cr,
 181                                           const gmx_reverse_top_t *rt,
 182                                           const char *moltypename,
 183                                           const reverse_ilist_t *ril,
 184                                           int a_start, int a_end,
 185                                           int nat_mol, int nmol,
 186                                           const t_idef *idef)
 187 {
 188     int *assigned;
 189     int  nril_mol = ril->index[nat_mol];
 190     snew(assigned, nmol*nril_mol);
 191
 192     int *gatindex = cr->dd->gatindex;
 193     for (int ftype = 0; ftype < F_NRE; ftype++)
 194     {
 195         if (dd_check_ftype(ftype, rt->bBCheck, rt->bConstr, rt->bSettle))
 196         {
 197             int            nral = NRAL(ftype);
 198             const t_ilist *il   = &idef->il[ftype];
 199             const t_iatom *ia   = il->iatoms;
 200             for (int i = 0; i < il->nr; i += 1+nral)
 201             {
 202                 int a0 = gatindex[ia[1]];
 203                 /* Check if this interaction is in
 204                  * the currently checked molblock.
 205                  */
 206                 if (a0 >= a_start && a0 < a_end)
 207                 {
 208                     int  mol    = (a0 - a_start)/nat_mol;
 209                     int  a0_mol = (a0 - a_start) - mol*nat_mol;
 210                     int  j_mol  = ril->index[a0_mol];
 211                     bool found  = false;
 212                     while (j_mol < ril->index[a0_mol+1] && !found)
 213                     {
 214                         int j       = mol*nril_mol + j_mol;
 215                         int ftype_j = ril->il[j_mol];
 216                         /* Here we need to check if this interaction has
 217                          * not already been assigned, since we could have
 218                          * multiply defined interactions.
 219                          */
 220                         if (ftype == ftype_j && ia[0] == ril->il[j_mol+1] &&
 221                             assigned[j] == 0)
 222                         {
 223                             /* Check the atoms */
 224                             found = true;
 225                             for (int a = 0; a < nral; a++)
 226                             {
 227                                 if (gatindex[ia[1+a]] !=
 228                                     a_start + mol*nat_mol + ril->il[j_mol+2+a])
 229                                 {
 230                                     found = false;
 231                                 }
 232                             }
 233                             if (found)
 234                             {
 235                                 assigned[j] = 1;
 236                             }
 237                         }
 238                         j_mol += 2 + nral_rt(ftype_j);
 239                     }
 240                     if (!found)
 241                     {
 242                         gmx_incons("Some interactions seem to be assigned multiple times");
 243                     }
 244                 }
 245                 ia += 1 + nral;
 246             }
 247         }
 248     }
 249
 250     gmx_sumi(nmol*nril_mol, assigned, cr);
 251
 252     int nprint = 10;
 253     int i      = 0;
 254     for (int mol = 0; mol < nmol; mol++)
 255     {
 256         int j_mol = 0;
 257         while (j_mol < nril_mol)
 258         {
 259             int ftype = ril->il[j_mol];
 260             int nral  = NRAL(ftype);
 261             int j     = mol*nril_mol + j_mol;
 262             if (assigned[j] == 0 &&
 263                 !(interaction_function[ftype].flags & IF_VSITE))
 264             {
 265                 if (DDMASTER(cr->dd))
 266                 {
 267                     if (i == 0)
 268                     {
 269                         print_error_header(fplog, moltypename, nprint);
 270                     }
 271                     fprintf(fplog, "%20s atoms",
 272                             interaction_function[ftype].longname);
 273                     fprintf(stderr, "%20s atoms",
 274                             interaction_function[ftype].longname);
 275                     int a;
 276                     for (a = 0; a < nral; a++)
 277                     {
 278                         fprintf(fplog, "%5d", ril->il[j_mol+2+a]+1);
 279                         fprintf(stderr, "%5d", ril->il[j_mol+2+a]+1);
 280                     }
 281                     while (a < 4)
 282                     {
 283                         fprintf(fplog, "     ");
 284                         fprintf(stderr, "     ");
 285                         a++;
 286                     }
 287                     fprintf(fplog, " global");
 288                     fprintf(stderr, " global");
 289                     for (a = 0; a < nral; a++)
 290                     {
 291                         fprintf(fplog, "%6d",
 292                                 a_start+mol*nat_mol+ril->il[j_mol+2+a]+1);
 293                         fprintf(stderr, "%6d",
 294                                 a_start+mol*nat_mol+ril->il[j_mol+2+a]+1);
 295                     }
 296                     fprintf(fplog, "\n");
 297                     fprintf(stderr, "\n");
 298                 }
 299                 i++;
 300                 if (i >= nprint)
 301                 {
 302                     break;
 303                 }
 304             }
 305             j_mol += 2 + nral_rt(ftype);
 306         }
 307     }
 308
 309     sfree(assigned);
 310 }
 311
 312 /*! \brief Help print error output when interactions are missing */
 313 static void print_missing_interactions_atoms(FILE *fplog, t_commrec *cr,
 314                                              const gmx_mtop_t *mtop,
 315                                              const t_idef *idef)
 316 {
 317     const gmx_reverse_top_t *rt = cr->dd->reverse_top;
 318
 319     /* Print the atoms in the missing interactions per molblock */
 320     int a_end = 0;
 321     for (const gmx_molblock_t &molb :  mtop->molblock)
 322     {
 323         int a_start = a_end;
 324         a_end       = a_start + molb.nmol*molb.natoms_mol;
 325
 326         print_missing_interactions_mb(fplog, cr, rt,
 327                                       *(mtop->moltype[molb.type].name),
 328                                       &rt->ril_mt[molb.type],
 329                                       a_start, a_end, molb.natoms_mol,
 330                                       molb.nmol,
 331                                       idef);
 332     }
 333 }
 334
 335 void dd_print_missing_interactions(FILE *fplog, t_commrec *cr,
 336                                    int local_count,
 337                                    const gmx_mtop_t *top_global,
 338                                    const gmx_localtop_t *top_local,
 339                                    t_state *state_local)
 340 {
 341     int             ndiff_tot, cl[F_NRE], n, ndiff, rest_global, rest_local;
 342     int             ftype, nral;
 343     char            buf[STRLEN];
 344     gmx_domdec_t   *dd;
 345
 346     dd = cr->dd;
 347
 348     if (fplog)
 349     {
 350         fprintf(fplog, "\nNot all bonded interactions have been properly assigned to the domain decomposition cells\n");
 351         fflush(fplog);
 352     }
 353
 354     ndiff_tot = local_count - dd->nbonded_global;
 355
 356     for (ftype = 0; ftype < F_NRE; ftype++)
 357     {
 358         nral      = NRAL(ftype);
 359         cl[ftype] = top_local->idef.il[ftype].nr/(1+nral);
 360     }
 361
 362     gmx_sumi(F_NRE, cl, cr);
 363
 364     if (DDMASTER(dd))
 365     {
 366         if (fplog)
 367         {
 368             fprintf(fplog, "\nA list of missing interactions:\n");
 369         }
 370         fprintf(stderr, "\nA list of missing interactions:\n");
 371         rest_global = dd->nbonded_global;
 372         rest_local  = local_count;
 373         for (ftype = 0; ftype < F_NRE; ftype++)
 374         {
 375             /* In the reverse and local top all constraints are merged
 376              * into F_CONSTR. So in the if statement we skip F_CONSTRNC
 377              * and add these constraints when doing F_CONSTR.
 378              */
 379             if (((interaction_function[ftype].flags & IF_BOND) &&
 380                  (dd->reverse_top->bBCheck
 381                   || !(interaction_function[ftype].flags & IF_LIMZERO)))
 382                 || (dd->reverse_top->bConstr && ftype == F_CONSTR)
 383                 || (dd->reverse_top->bSettle && ftype == F_SETTLE))
 384             {
 385                 n    = gmx_mtop_ftype_count(top_global, ftype);
 386                 if (ftype == F_CONSTR)
 387                 {
 388                     n += gmx_mtop_ftype_count(top_global, F_CONSTRNC);
 389                 }
 390                 ndiff = cl[ftype] - n;
 391                 if (ndiff != 0)
 392                 {
 393                     sprintf(buf, "%20s of %6d missing %6d",
 394                             interaction_function[ftype].longname, n, -ndiff);
 395                     if (fplog)
 396                     {
 397                         fprintf(fplog, "%s\n", buf);
 398                     }
 399                     fprintf(stderr, "%s\n", buf);
 400                 }
 401                 rest_global -= n;
 402                 rest_local  -= cl[ftype];
 403             }
 404         }
 405
 406         ndiff = rest_local - rest_global;
 407         if (ndiff != 0)
 408         {
 409             sprintf(buf, "%20s of %6d missing %6d", "exclusions",
 410                     rest_global, -ndiff);
 411             if (fplog)
 412             {
 413                 fprintf(fplog, "%s\n", buf);
 414             }
 415             fprintf(stderr, "%s\n", buf);
 416         }
 417     }
 418
 419     print_missing_interactions_atoms(fplog, cr, top_global, &top_local->idef);
 420     write_dd_pdb("dd_dump_err", 0, "dump", top_global, cr,
 421                  -1, as_rvec_array(state_local->x.data()), state_local->box);
 422
 423     std::string errorMessage;
 424
 425     if (ndiff_tot > 0)
 426     {
 427         errorMessage = "One or more interactions were assigned to multiple domains of the domain decompostion. Please report this bug.";
 428     }
 429     else
 430     {
 431         errorMessage = gmx::formatString("%d of the %d bonded interactions could not be calculated because some atoms involved moved further apart than the multi-body cut-off distance (%g nm) or the two-body cut-off distance (%g nm), see option -rdd, for pairs and tabulated bonds also see option -ddcheck", -ndiff_tot, cr->dd->nbonded_global, dd_cutoff_multibody(dd), dd_cutoff_twobody(dd));
 432     }
 433     gmx_fatal_collective(FARGS, cr->mpi_comm_mygroup, MASTER(cr), errorMessage.c_str());
 434 }
 435
 436 /*! \brief Return global topology molecule information for global atom index \p i_gl */
 437 static void global_atomnr_to_moltype_ind(gmx_reverse_top_t *rt, int i_gl,
 438                                          int *mb, int *mt, int *mol, int *i_mol)
 439 {
 440     molblock_ind_t *mbi   = rt->mbi;
 441     int             start = 0;
 442     int             end   = rt->nmolblock; /* exclusive */
 443     int             mid;
 444
 445     /* binary search for molblock_ind */
 446     while (TRUE)
 447     {
 448         mid = (start+end)/2;
 449         if (i_gl >= mbi[mid].a_end)
 450         {
 451             start = mid+1;
 452         }
 453         else if (i_gl < mbi[mid].a_start)
 454         {
 455             end = mid;
 456         }
 457         else
 458         {
 459             break;
 460         }
 461     }
 462
 463     *mb  = mid;
 464     mbi += mid;
 465
 466     *mt    = mbi->type;
 467     *mol   = (i_gl - mbi->a_start) / mbi->natoms_mol;
 468     *i_mol = (i_gl - mbi->a_start) - (*mol)*mbi->natoms_mol;
 469 }
 470
 471 /*! \brief Count the exclusions for all atoms in \p cgs */
 472 static void count_excls(const t_block *cgs, const t_blocka *excls,
 473                         int *n_excl, int *n_intercg_excl, int *n_excl_at_max)
 474 {
 475     int cg, at0, at1, at, excl, atj;
 476
 477     *n_excl         = 0;
 478     *n_intercg_excl = 0;
 479     *n_excl_at_max  = 0;
 480     for (cg = 0; cg < cgs->nr; cg++)
 481     {
 482         at0 = cgs->index[cg];
 483         at1 = cgs->index[cg+1];
 484         for (at = at0; at < at1; at++)
 485         {
 486             for (excl = excls->index[at]; excl < excls->index[at+1]; excl++)
 487             {
 488                 atj = excls->a[excl];
 489                 if (atj > at)
 490                 {
 491                     (*n_excl)++;
 492                     if (atj < at0 || atj >= at1)
 493                     {
 494                         (*n_intercg_excl)++;
 495                     }
 496                 }
 497             }
 498
 499             *n_excl_at_max = std::max(*n_excl_at_max,
 500                                       excls->index[at+1] - excls->index[at]);
 501         }
 502     }
 503 }
 504
 505 /*! \brief Run the reverse ilist generation and store it in r_il when \p bAssign = TRUE */
 506 static int low_make_reverse_ilist(const t_ilist *il_mt, const t_atom *atom,
 507                                   const int * const * vsite_pbc,
 508                                   int *count,
 509                                   gmx_bool bConstr, gmx_bool bSettle,
 510                                   gmx_bool bBCheck,
 511                                   int *r_index, int *r_il,
 512                                   gmx_bool bLinkToAllAtoms,
 513                                   gmx_bool bAssign)
 514 {
 515     int            ftype, nral, i, j, nlink, link;
 516     const t_ilist *il;
 517     const t_iatom *ia;
 518     int            a;
 519     int            nint;
 520     gmx_bool       bVSite;
 521
 522     nint = 0;
 523     for (ftype = 0; ftype < F_NRE; ftype++)
 524     {
 525         if ((interaction_function[ftype].flags & (IF_BOND | IF_VSITE)) ||
 526             (bConstr && (ftype == F_CONSTR || ftype == F_CONSTRNC)) ||
 527             (bSettle && ftype == F_SETTLE))
 528         {
 529             bVSite = (interaction_function[ftype].flags & IF_VSITE);
 530             nral   = NRAL(ftype);
 531             il     = &il_mt[ftype];
 532             for (i = 0; i < il->nr; i += 1+nral)
 533             {
 534                 ia = il->iatoms + i;
 535                 if (bLinkToAllAtoms)
 536                 {
 537                     if (bVSite)
 538                     {
 539                         /* We don't need the virtual sites for the cg-links */
 540                         nlink = 0;
 541                     }
 542                     else
 543                     {
 544                         nlink = nral;
 545                     }
 546                 }
 547                 else
 548                 {
 549                     /* Couple to the first atom in the interaction */
 550                     nlink = 1;
 551                 }
 552                 for (link = 0; link < nlink; link++)
 553                 {
 554                     a = ia[1+link];
 555                     if (bAssign)
 556                     {
 557                         assert(r_il); //with bAssign not allowed to be null
 558                         assert(r_index);
 559                         r_il[r_index[a]+count[a]] =
 560                             (ftype == F_CONSTRNC ? F_CONSTR : ftype);
 561                         r_il[r_index[a]+count[a]+1] = ia[0];
 562                         for (j = 1; j < 1+nral; j++)
 563                         {
 564                             /* Store the molecular atom number */
 565                             r_il[r_index[a]+count[a]+1+j] = ia[j];
 566                         }
 567                     }
 568                     if (interaction_function[ftype].flags & IF_VSITE)
 569                     {
 570                         if (bAssign)
 571                         {
 572                             /* Add an entry to iatoms for storing
 573                              * which of the constructing atoms are
 574                              * vsites again.
 575                              */
 576                             r_il[r_index[a]+count[a]+2+nral] = 0;
 577                             for (j = 2; j < 1+nral; j++)
 578                             {
 579                                 if (atom[ia[j]].ptype == eptVSite)
 580                                 {
 581                                     r_il[r_index[a]+count[a]+2+nral] |= (2<<j);
 582                                 }
 583                             }
 584                             /* Store vsite pbc atom in a second extra entry */
 585                             r_il[r_index[a]+count[a]+2+nral+1] =
 586                                 (vsite_pbc ? vsite_pbc[ftype-F_VSITE2][i/(1+nral)] : -2);
 587                         }
 588                     }
 589                     else
 590                     {
 591                         /* We do not count vsites since they are always
 592                          * uniquely assigned and can be assigned
 593                          * to multiple nodes with recursive vsites.
 594                          */
 595                         if (bBCheck ||
 596                             !(interaction_function[ftype].flags & IF_LIMZERO))
 597                         {
 598                             nint++;
 599                         }
 600                     }
 601                     count[a] += 2 + nral_rt(ftype);
 602                 }
 603             }
 604         }
 605     }
 606
 607     return nint;
 608 }
 609
 610 /*! \brief Make the reverse ilist: a list of bonded interactions linked to atoms */
 611 static int make_reverse_ilist(const t_ilist *ilist,
 612                               const t_atoms *atoms,
 613                               const int * const * vsite_pbc,
 614                               gmx_bool bConstr, gmx_bool bSettle,
 615                               gmx_bool bBCheck,
 616                               gmx_bool bLinkToAllAtoms,
 617                               reverse_ilist_t *ril_mt)
 618 {
 619     int nat_mt, *count, i, nint_mt;
 620
 621     /* Count the interactions */
 622     nat_mt = atoms->nr;
 623     snew(count, nat_mt);
 624     low_make_reverse_ilist(ilist, atoms->atom, vsite_pbc,
 625                            count,
 626                            bConstr, bSettle, bBCheck, nullptr, nullptr,
 627                            bLinkToAllAtoms, FALSE);
 628
 629     snew(ril_mt->index, nat_mt+1);
 630     ril_mt->index[0] = 0;
 631     for (i = 0; i < nat_mt; i++)
 632     {
 633         ril_mt->index[i+1] = ril_mt->index[i] + count[i];
 634         count[i]           = 0;
 635     }
 636     snew(ril_mt->il, ril_mt->index[nat_mt]);
 637
 638     /* Store the interactions */
 639     nint_mt =
 640         low_make_reverse_ilist(ilist, atoms->atom, vsite_pbc,
 641                                count,
 642                                bConstr, bSettle, bBCheck,
 643                                ril_mt->index, ril_mt->il,
 644                                bLinkToAllAtoms, TRUE);
 645
 646     sfree(count);
 647
 648     return nint_mt;
 649 }
 650
 651 /*! \brief Destroys a reverse_ilist_t struct */
 652 static void destroy_reverse_ilist(reverse_ilist_t *ril)
 653 {
 654     sfree(ril->index);
 655     sfree(ril->il);
 656 }
 657
 658 /*! \brief Generate the reverse topology */
 659 static gmx_reverse_top_t *make_reverse_top(const gmx_mtop_t *mtop, gmx_bool bFE,
 660                                            const int * const * const * vsite_pbc_molt,
 661                                            gmx_bool bConstr, gmx_bool bSettle,
 662                                            gmx_bool bBCheck, int *nint)
 663 {
 664     gmx_reverse_top_t *rt;
 665     int               *nint_mt;
 666     int                thread;
 667
 668     snew(rt, 1);
 669
 670     /* Should we include constraints (for SHAKE) in rt? */
 671     rt->bConstr = bConstr;
 672     rt->bSettle = bSettle;
 673     rt->bBCheck = bBCheck;
 674
 675     rt->bInterCGInteractions = mtop->bIntermolecularInteractions;
 676     snew(nint_mt, mtop->moltype.size());
 677     snew(rt->ril_mt, mtop->moltype.size());
 678     rt->ril_mt_tot_size = 0;
 679     for (size_t mt = 0; mt < mtop->moltype.size(); mt++)
 680     {
 681         const gmx_moltype_t &molt = mtop->moltype[mt];
 682         if (molt.cgs.nr > 1)
 683         {
 684             rt->bInterCGInteractions = TRUE;
 685         }
 686
 687         /* Make the atom to interaction list for this molecule type */
 688         nint_mt[mt] =
 689             make_reverse_ilist(molt.ilist, &molt.atoms,
 690                                vsite_pbc_molt ? vsite_pbc_molt[mt] : nullptr,
 691                                rt->bConstr, rt->bSettle, rt->bBCheck, FALSE,
 692                                &rt->ril_mt[mt]);
 693
 694         rt->ril_mt_tot_size += rt->ril_mt[mt].index[molt.atoms.nr];
 695     }
 696     if (debug)
 697     {
 698         fprintf(debug, "The total size of the atom to interaction index is %d integers\n", rt->ril_mt_tot_size);
 699     }
 700
 701     *nint = 0;
 702     for (const gmx_molblock_t &molblock : mtop->molblock)
 703     {
 704         *nint += molblock.nmol*nint_mt[molblock.type];
 705     }
 706     sfree(nint_mt);
 707
 708     /* Make an intermolecular reverse top, if necessary */
 709     rt->bIntermolecularInteractions = mtop->bIntermolecularInteractions;
 710     if (rt->bIntermolecularInteractions)
 711     {
 712         t_atoms atoms_global;
 713
 714         rt->ril_intermol.index = nullptr;
 715         rt->ril_intermol.il    = nullptr;
 716
 717         atoms_global.nr   = mtop->natoms;
 718         atoms_global.atom = nullptr; /* Only used with virtual sites */
 719
 720         *nint +=
 721             make_reverse_ilist(mtop->intermolecular_ilist, &atoms_global,
 722                                nullptr,
 723                                rt->bConstr, rt->bSettle, rt->bBCheck, FALSE,
 724                                &rt->ril_intermol);
 725     }
 726
 727     if (bFE && gmx_mtop_bondeds_free_energy(mtop))
 728     {
 729         rt->ilsort = ilsortFE_UNSORTED;
 730     }
 731     else
 732     {
 733         rt->ilsort = ilsortNO_FE;
 734     }
 735
 736     /* Make a molblock index for fast searching */
 737     snew(rt->mbi, mtop->molblock.size());
 738     rt->nmolblock = mtop->molblock.size();
 739     int i         = 0;
 740     for (size_t mb = 0; mb < mtop->molblock.size(); mb++)
 741     {
 742         rt->mbi[mb].a_start    = i;
 743         i += mtop->molblock[mb].nmol*mtop->molblock[mb].natoms_mol;
 744         rt->mbi[mb].a_end      = i;
 745         rt->mbi[mb].natoms_mol = mtop->molblock[mb].natoms_mol;
 746         rt->mbi[mb].type       = mtop->molblock[mb].type;
 747     }
 748
 749     rt->nthread = gmx_omp_nthreads_get(emntDomdec);
 750     snew(rt->th_work, rt->nthread);
 751     if (vsite_pbc_molt != nullptr)
 752     {
 753         for (thread = 0; thread < rt->nthread; thread++)
 754         {
 755             snew(rt->th_work[thread].vsite_pbc, F_VSITEN-F_VSITE2+1);
 756             snew(rt->th_work[thread].vsite_pbc_nalloc, F_VSITEN-F_VSITE2+1);
 757         }
 758     }
 759
 760     return rt;
 761 }
 762
 763 void dd_make_reverse_top(FILE *fplog,
 764                          gmx_domdec_t *dd, const gmx_mtop_t *mtop,
 765                          const gmx_vsite_t *vsite,
 766                          const t_inputrec *ir, gmx_bool bBCheck)
 767 {
 768     if (fplog)
 769     {
 770         fprintf(fplog, "\nLinking all bonded interactions to atoms\n");
 771     }
 772
 773     /* If normal and/or settle constraints act only within charge groups,
 774      * we can store them in the reverse top and simply assign them to domains.
 775      * Otherwise we need to assign them to multiple domains and set up
 776      * the parallel version constraint algorithm(s).
 777      */
 778
 779     dd->reverse_top = make_reverse_top(mtop, ir->efep != efepNO,
 780                                        vsite ? vsite->vsite_pbc_molt : nullptr,
 781                                        !dd->bInterCGcons, !dd->bInterCGsettles,
 782                                        bBCheck, &dd->nbonded_global);
 783
 784     gmx_reverse_top_t *rt = dd->reverse_top;
 785
 786     /* With the Verlet scheme, exclusions are handled in the non-bonded
 787      * kernels and only exclusions inside the cut-off lead to exclusion
 788      * forces. Since each atom pair is treated at most once in the non-bonded
 789      * kernels, it doesn't matter if the exclusions for the same atom pair
 790      * appear multiple times in the exclusion list. In contrast, the, old,
 791      * group cut-off scheme loops over a list of exclusions, so there each
 792      * excluded pair should appear exactly once.
 793      */
 794     rt->bExclRequired = (ir->cutoff_scheme == ecutsGROUP &&
 795                          inputrecExclForces(ir));
 796
 797     int nexcl          = 0;
 798     dd->n_intercg_excl = 0;
 799     rt->n_excl_at_max  = 0;
 800     for (const gmx_molblock_t &molb : mtop->molblock)
 801     {
 802         int                  n_excl_mol, n_excl_icg, n_excl_at_max;
 803
 804         const gmx_moltype_t &molt = mtop->moltype[molb.type];
 805         count_excls(&molt.cgs, &molt.excls,
 806                     &n_excl_mol, &n_excl_icg, &n_excl_at_max);
 807         nexcl              += molb.nmol*n_excl_mol;
 808         dd->n_intercg_excl += molb.nmol*n_excl_icg;
 809         rt->n_excl_at_max   = std::max(rt->n_excl_at_max, n_excl_at_max);
 810     }
 811     if (rt->bExclRequired)
 812     {
 813         dd->nbonded_global += nexcl;
 814         if (EEL_FULL(ir->coulombtype) && dd->n_intercg_excl > 0 && fplog)
 815         {
 816             fprintf(fplog, "There are %d inter charge-group exclusions,\n"
 817                     "will use an extra communication step for exclusion forces for %s\n",
 818                     dd->n_intercg_excl, eel_names[ir->coulombtype]);
 819         }
 820     }
 821
 822     if (vsite && vsite->n_intercg_vsite > 0)
 823     {
 824         if (fplog)
 825         {
 826             fprintf(fplog, "There are %d inter charge-group virtual sites,\n"
 827                     "will an extra communication step for selected coordinates and forces\n",
 828                     vsite->n_intercg_vsite);
 829         }
 830         init_domdec_vsites(dd, vsite->n_intercg_vsite);
 831     }
 832
 833     if (dd->bInterCGcons || dd->bInterCGsettles)
 834     {
 835         init_domdec_constraints(dd, mtop);
 836     }
 837     if (fplog)
 838     {
 839         fprintf(fplog, "\n");
 840     }
 841 }
 842
 843 /*! \brief Store a vsite interaction at the end of \p il
 844  *
 845  * This routine is very similar to add_ifunc, but vsites interactions
 846  * have more work to do than other kinds of interactions, and the
 847  * complex way nral (and thus vector contents) depends on ftype
 848  * confuses static analysis tools unless we fuse the vsite
 849  * atom-indexing organization code with the ifunc-adding code, so that
 850  * they can see that nral is the same value. */
 851 static inline void
 852 add_ifunc_for_vsites(t_iatom *tiatoms, gmx_ga2la_t *ga2la,
 853                      int nral, gmx_bool bHomeA,
 854                      int a, int a_gl, int a_mol,
 855                      const t_iatom *iatoms,
 856                      t_ilist *il)
 857 {
 858     t_iatom *liatoms;
 859
 860     if (il->nr+1+nral > il->nalloc)
 861     {
 862         il->nalloc = over_alloc_large(il->nr+1+nral);
 863         srenew(il->iatoms, il->nalloc);
 864     }
 865     liatoms = il->iatoms + il->nr;
 866     il->nr += 1 + nral;
 867
 868     /* Copy the type */
 869     tiatoms[0] = iatoms[0];
 870
 871     if (bHomeA)
 872     {
 873         /* We know the local index of the first atom */
 874         tiatoms[1] = a;
 875     }
 876     else
 877     {
 878         /* Convert later in make_local_vsites */
 879         tiatoms[1] = -a_gl - 1;
 880     }
 881
 882     for (int k = 2; k < 1+nral; k++)
 883     {
 884         int ak_gl = a_gl + iatoms[k] - a_mol;
 885         if (!ga2la_get_home(ga2la, ak_gl, &tiatoms[k]))
 886         {
 887             /* Copy the global index, convert later in make_local_vsites */
 888             tiatoms[k] = -(ak_gl + 1);
 889         }
 890         // Note that ga2la_get_home always sets the third parameter if
 891         // it returns TRUE
 892     }
 893     for (int k = 0; k < 1+nral; k++)
 894     {
 895         liatoms[k] = tiatoms[k];
 896     }
 897 }
 898
 899 /*! \brief Store a bonded interaction at the end of \p il */
 900 static inline void add_ifunc(int nral, t_iatom *tiatoms, t_ilist *il)
 901 {
 902     t_iatom *liatoms;
 903     int      k;
 904
 905     if (il->nr+1+nral > il->nalloc)
 906     {
 907         il->nalloc = over_alloc_large(il->nr+1+nral);
 908         srenew(il->iatoms, il->nalloc);
 909     }
 910     liatoms = il->iatoms + il->nr;
 911     for (k = 0; k <= nral; k++)
 912     {
 913         liatoms[k] = tiatoms[k];
 914     }
 915     il->nr += 1 + nral;
 916 }
 917
 918 /*! \brief Store a position restraint in idef and iatoms, complex because the parameters are different for each entry */
 919 static void add_posres(int mol, int a_mol, const gmx_molblock_t *molb,
 920                        t_iatom *iatoms, const t_iparams *ip_in,
 921                        t_idef *idef)
 922 {
 923     int        n, a_molb;
 924     t_iparams *ip;
 925
 926     /* This position restraint has not been added yet,
 927      * so it's index is the current number of position restraints.
 928      */
 929     n = idef->il[F_POSRES].nr/2;
 930     if (n+1 > idef->iparams_posres_nalloc)
 931     {
 932         idef->iparams_posres_nalloc = over_alloc_dd(n+1);
 933         srenew(idef->iparams_posres, idef->iparams_posres_nalloc);
 934     }
 935     ip = &idef->iparams_posres[n];
 936     /* Copy the force constants */
 937     *ip = ip_in[iatoms[0]];
 938
 939     /* Get the position restraint coordinates from the molblock */
 940     a_molb = mol*molb->natoms_mol + a_mol;
 941     if (a_molb >= molb->nposres_xA)
 942     {
 943         gmx_incons("Not enough position restraint coordinates");
 944     }
 945     ip->posres.pos0A[XX] = molb->posres_xA[a_molb][XX];
 946     ip->posres.pos0A[YY] = molb->posres_xA[a_molb][YY];
 947     ip->posres.pos0A[ZZ] = molb->posres_xA[a_molb][ZZ];
 948     if (molb->nposres_xB > 0)
 949     {
 950         ip->posres.pos0B[XX] = molb->posres_xB[a_molb][XX];
 951         ip->posres.pos0B[YY] = molb->posres_xB[a_molb][YY];
 952         ip->posres.pos0B[ZZ] = molb->posres_xB[a_molb][ZZ];
 953     }
 954     else
 955     {
 956         ip->posres.pos0B[XX] = ip->posres.pos0A[XX];
 957         ip->posres.pos0B[YY] = ip->posres.pos0A[YY];
 958         ip->posres.pos0B[ZZ] = ip->posres.pos0A[ZZ];
 959     }
 960     /* Set the parameter index for idef->iparams_posre */
 961     iatoms[0] = n;
 962 }
 963
 964 /*! \brief Store a flat-bottomed position restraint in idef and iatoms, complex because the parameters are different for each entry */
 965 static void add_fbposres(int mol, int a_mol, const gmx_molblock_t *molb,
 966                          t_iatom *iatoms, const t_iparams *ip_in,
 967                          t_idef *idef)
 968 {
 969     int        n, a_molb;
 970     t_iparams *ip;
 971
 972     /* This flat-bottom position restraint has not been added yet,
 973      * so it's index is the current number of position restraints.
 974      */
 975     n = idef->il[F_FBPOSRES].nr/2;
 976     if (n+1 > idef->iparams_fbposres_nalloc)
 977     {
 978         idef->iparams_fbposres_nalloc = over_alloc_dd(n+1);
 979         srenew(idef->iparams_fbposres, idef->iparams_fbposres_nalloc);
 980     }
 981     ip = &idef->iparams_fbposres[n];
 982     /* Copy the force constants */
 983     *ip = ip_in[iatoms[0]];
 984
 985     /* Get the position restriant coordinats from the molblock */
 986     a_molb = mol*molb->natoms_mol + a_mol;
 987     if (a_molb >= molb->nposres_xA)
 988     {
 989         gmx_incons("Not enough position restraint coordinates");
 990     }
 991     /* Take reference positions from A position of normal posres */
 992     ip->fbposres.pos0[XX] = molb->posres_xA[a_molb][XX];
 993     ip->fbposres.pos0[YY] = molb->posres_xA[a_molb][YY];
 994     ip->fbposres.pos0[ZZ] = molb->posres_xA[a_molb][ZZ];
 995
 996     /* Note: no B-type for flat-bottom posres */
 997
 998     /* Set the parameter index for idef->iparams_posre */
 999     iatoms[0] = n;
1000 }
1001
1002 /*! \brief Store a virtual site interaction, complex because of PBC and recursion */
1003 static void add_vsite(gmx_ga2la_t *ga2la, const int *index, const int *rtil,
1004                       int ftype, int nral,
1005                       gmx_bool bHomeA, int a, int a_gl, int a_mol,
1006                       const t_iatom *iatoms,
1007                       t_idef *idef, int **vsite_pbc, int *vsite_pbc_nalloc)
1008 {
1009     int     k, vsi, pbc_a_mol;
1010     t_iatom tiatoms[1+MAXATOMLIST];
1011     int     j, ftype_r, nral_r;
1012
1013     /* Add this interaction to the local topology */
1014     add_ifunc_for_vsites(tiatoms, ga2la, nral, bHomeA, a, a_gl, a_mol, iatoms, &idef->il[ftype]);
1015
1016     if (vsite_pbc)
1017     {
1018         vsi = idef->il[ftype].nr/(1+nral) - 1;
1019         if (vsi >= vsite_pbc_nalloc[ftype-F_VSITE2])
1020         {
1021             vsite_pbc_nalloc[ftype-F_VSITE2] = over_alloc_large(vsi+1);
1022             srenew(vsite_pbc[ftype-F_VSITE2], vsite_pbc_nalloc[ftype-F_VSITE2]);
1023         }
1024         if (bHomeA)
1025         {
1026             pbc_a_mol = iatoms[1+nral+1];
1027             if (pbc_a_mol < 0)
1028             {
1029                 /* The pbc flag is one of the following two options:
1030                  * -2: vsite and all constructing atoms are within the same cg, no pbc
1031                  * -1: vsite and its first constructing atom are in the same cg, do pbc
1032                  */
1033                 vsite_pbc[ftype-F_VSITE2][vsi] = pbc_a_mol;
1034             }
1035             else
1036             {
1037                 /* Set the pbc atom for this vsite so we can make its pbc
1038                  * identical to the rest of the atoms in its charge group.
1039                  * Since the order of the atoms does not change within a charge
1040                  * group, we do not need the global to local atom index.
1041                  */
1042                 vsite_pbc[ftype-F_VSITE2][vsi] = a + pbc_a_mol - iatoms[1];
1043             }
1044         }
1045         else
1046         {
1047             /* This vsite is non-home (required for recursion),
1048              * and therefore there is no charge group to match pbc with.
1049              * But we always turn on full_pbc to assure that higher order
1050              * recursion works correctly.
1051              */
1052             vsite_pbc[ftype-F_VSITE2][vsi] = -1;
1053         }
1054     }
1055
1056     if (iatoms[1+nral])
1057     {
1058         /* Check for recursion */
1059         for (k = 2; k < 1+nral; k++)
1060         {
1061             if ((iatoms[1+nral] & (2<<k)) && (tiatoms[k] < 0))
1062             {
1063                 /* This construction atoms is a vsite and not a home atom */
1064                 if (gmx_debug_at)
1065                 {
1066                     fprintf(debug, "Constructing atom %d of vsite atom %d is a vsite and non-home\n", iatoms[k]+1, a_mol+1);
1067                 }
1068                 /* Find the vsite construction */
1069
1070                 /* Check all interactions assigned to this atom */
1071                 j = index[iatoms[k]];
1072                 while (j < index[iatoms[k]+1])
1073                 {
1074                     ftype_r = rtil[j++];
1075                     nral_r  = NRAL(ftype_r);
1076                     if (interaction_function[ftype_r].flags & IF_VSITE)
1077                     {
1078                         /* Add this vsite (recursion) */
1079                         add_vsite(ga2la, index, rtil, ftype_r, nral_r,
1080                                   FALSE, -1, a_gl+iatoms[k]-iatoms[1], iatoms[k],
1081                                   rtil+j, idef, vsite_pbc, vsite_pbc_nalloc);
1082                         j += 1 + nral_r + 2;
1083                     }
1084                     else
1085                     {
1086                         j += 1 + nral_r;
1087                     }
1088                 }
1089             }
1090         }
1091     }
1092 }
1093
1094 /*! \brief Update the local atom to local charge group index */
1095 static void make_la2lc(gmx_domdec_t *dd)
1096 {
1097     int *cgindex, *la2lc, cg, a;
1098
1099     cgindex = dd->cgindex;
1100
1101     if (dd->nat_tot > dd->la2lc_nalloc)
1102     {
1103         dd->la2lc_nalloc = over_alloc_dd(dd->nat_tot);
1104         snew(dd->la2lc, dd->la2lc_nalloc);
1105     }
1106     la2lc = dd->la2lc;
1107
1108     /* Make the local atom to local cg index */
1109     for (cg = 0; cg < dd->ncg_tot; cg++)
1110     {
1111         for (a = cgindex[cg]; a < cgindex[cg+1]; a++)
1112         {
1113             la2lc[a] = cg;
1114         }
1115     }
1116 }
1117
1118 /*! \brief Returns the squared distance between charge groups \p i and \p j */
1119 static real dd_dist2(t_pbc *pbc_null, rvec *cg_cm, const int *la2lc, int i, int j)
1120 {
1121     rvec dx;
1122
1123     if (pbc_null)
1124     {
1125         pbc_dx_aiuc(pbc_null, cg_cm[la2lc[i]], cg_cm[la2lc[j]], dx);
1126     }
1127     else
1128     {
1129         rvec_sub(cg_cm[la2lc[i]], cg_cm[la2lc[j]], dx);
1130     }
1131
1132     return norm2(dx);
1133 }
1134
1135 /*! \brief Append t_blocka block structures 1 to nsrc in src to *dest */
1136 static void combine_blocka(t_blocka *dest, const thread_work_t *src, int nsrc)
1137 {
1138     int ni, na, s, i;
1139
1140     ni = src[nsrc-1].excl.nr;
1141     na = 0;
1142     for (s = 0; s < nsrc; s++)
1143     {
1144         na += src[s].excl.nra;
1145     }
1146     if (ni + 1 > dest->nalloc_index)
1147     {
1148         dest->nalloc_index = over_alloc_large(ni+1);
1149         srenew(dest->index, dest->nalloc_index);
1150     }
1151     if (dest->nra + na > dest->nalloc_a)
1152     {
1153         dest->nalloc_a = over_alloc_large(dest->nra+na);
1154         srenew(dest->a, dest->nalloc_a);
1155     }
1156     for (s = 1; s < nsrc; s++)
1157     {
1158         for (i = dest->nr+1; i < src[s].excl.nr+1; i++)
1159         {
1160             dest->index[i] = dest->nra + src[s].excl.index[i];
1161         }
1162         for (i = 0; i < src[s].excl.nra; i++)
1163         {
1164             dest->a[dest->nra+i] = src[s].excl.a[i];
1165         }
1166         dest->nr   = src[s].excl.nr;
1167         dest->nra += src[s].excl.nra;
1168     }
1169 }
1170
1171 /*! \brief Append t_idef structures 1 to nsrc in src to *dest,
1172  * virtual sites need special attention, as pbc info differs per vsite.
1173  */
1174 static void combine_idef(t_idef *dest, const thread_work_t *src, int nsrc,
1175                          gmx_vsite_t *vsite)
1176 {
1177     int ftype;
1178
1179     for (ftype = 0; ftype < F_NRE; ftype++)
1180     {
1181         int n, s;
1182
1183         n = 0;
1184         for (s = 1; s < nsrc; s++)
1185         {
1186             n += src[s].idef.il[ftype].nr;
1187         }
1188         if (n > 0)
1189         {
1190             t_ilist *ild;
1191
1192             ild = &dest->il[ftype];
1193
1194             if (ild->nr + n > ild->nalloc)
1195             {
1196                 ild->nalloc = over_alloc_large(ild->nr+n);
1197                 srenew(ild->iatoms, ild->nalloc);
1198             }
1199
1200             gmx_bool vpbc;
1201             int      nral1 = 0, ftv = 0;
1202
1203             vpbc = ((interaction_function[ftype].flags & IF_VSITE) &&
1204                     vsite->vsite_pbc_loc != nullptr);
1205             if (vpbc)
1206             {
1207                 nral1 = 1 + NRAL(ftype);
1208                 ftv   = ftype - F_VSITE2;
1209                 if ((ild->nr + n)/nral1 > vsite->vsite_pbc_loc_nalloc[ftv])
1210                 {
1211                     vsite->vsite_pbc_loc_nalloc[ftv] =
1212                         over_alloc_large((ild->nr + n)/nral1);
1213                     srenew(vsite->vsite_pbc_loc[ftv],
1214                            vsite->vsite_pbc_loc_nalloc[ftv]);
1215                 }
1216             }
1217
1218             for (s = 1; s < nsrc; s++)
1219             {
1220                 const t_ilist *ils;
1221                 int            i;
1222
1223                 ils = &src[s].idef.il[ftype];
1224                 for (i = 0; i < ils->nr; i++)
1225                 {
1226                     ild->iatoms[ild->nr+i] = ils->iatoms[i];
1227                 }
1228                 if (vpbc)
1229                 {
1230                     for (i = 0; i < ils->nr; i += nral1)
1231                     {
1232                         vsite->vsite_pbc_loc[ftv][(ild->nr+i)/nral1] =
1233                             src[s].vsite_pbc[ftv][i/nral1];
1234                     }
1235                 }
1236
1237                 ild->nr += ils->nr;
1238             }
1239
1240             /* Position restraints need an additional treatment */
1241             if (ftype == F_POSRES || ftype == F_FBPOSRES)
1242             {
1243                 int          nposres       = dest->il[ftype].nr/2;
1244                 // TODO: Simplify this code using std::vector
1245                 t_iparams * &iparams_dest  = (ftype == F_POSRES ? dest->iparams_posres : dest->iparams_fbposres);
1246                 int         &posres_nalloc = (ftype == F_POSRES ? dest->iparams_posres_nalloc : dest->iparams_fbposres_nalloc);
1247                 if (nposres > posres_nalloc)
1248                 {
1249                     posres_nalloc = over_alloc_large(nposres);
1250                     srenew(iparams_dest, posres_nalloc);
1251                 }
1252
1253                 /* Set nposres to the number of original position restraints in dest */
1254                 for (int s = 1; s < nsrc; s++)
1255                 {
1256                     nposres -= src[s].idef.il[ftype].nr/2;
1257                 }
1258
1259                 for (int s = 1; s < nsrc; s++)
1260                 {
1261                     const t_iparams *iparams_src = (ftype == F_POSRES ? src[s].idef.iparams_posres : src[s].idef.iparams_fbposres);
1262
1263                     for (int i = 0; i < src[s].idef.il[ftype].nr/2; i++)
1264                     {
1265                         /* Correct the index into iparams_posres */
1266                         dest->il[ftype].iatoms[nposres*2] = nposres;
1267                         /* Copy the position restraint force parameters */
1268                         iparams_dest[nposres]             = iparams_src[i];
1269                         nposres++;
1270                     }
1271                 }
1272             }
1273         }
1274     }
1275 }
1276
1277 /*! \brief Check and when available assign bonded interactions for local atom i
1278  */
1279 static inline void
1280 check_assign_interactions_atom(int i, int i_gl,
1281                                int mol, int i_mol,
1282                                const int *index, const int *rtil,
1283                                gmx_bool bInterMolInteractions,
1284                                int ind_start, int ind_end,
1285                                const gmx_domdec_t *dd,
1286                                const gmx_domdec_zones_t *zones,
1287                                const gmx_molblock_t *molb,
1288                                gmx_bool bRCheckMB, ivec rcheck, gmx_bool bRCheck2B,
1289                                real rc2,
1290                                int *la2lc,
1291                                t_pbc *pbc_null,
1292                                rvec *cg_cm,
1293                                const t_iparams *ip_in,
1294                                t_idef *idef,
1295                                int **vsite_pbc, int *vsite_pbc_nalloc,
1296                                int iz,
1297                                gmx_bool bBCheck,
1298                                int *nbonded_local)
1299 {
1300     int j;
1301
1302     j = ind_start;
1303     while (j < ind_end)
1304     {
1305         int            ftype;
1306         const t_iatom *iatoms;
1307         int            nral;
1308         t_iatom        tiatoms[1 + MAXATOMLIST];
1309
1310         ftype  = rtil[j++];
1311         iatoms = rtil + j;
1312         nral   = NRAL(ftype);
1313         if (ftype == F_SETTLE)
1314         {
1315             /* Settles are only in the reverse top when they
1316              * operate within a charge group. So we can assign
1317              * them without checks. We do this only for performance
1318              * reasons; it could be handled by the code below.
1319              */
1320             if (iz == 0)
1321             {
1322                 /* Home zone: add this settle to the local topology */
1323                 tiatoms[0] = iatoms[0];
1324                 tiatoms[1] = i;
1325                 tiatoms[2] = i + iatoms[2] - iatoms[1];
1326                 tiatoms[3] = i + iatoms[3] - iatoms[1];
1327                 add_ifunc(nral, tiatoms, &idef->il[ftype]);
1328                 (*nbonded_local)++;
1329             }
1330             j += 1 + nral;
1331         }
1332         else if (interaction_function[ftype].flags & IF_VSITE)
1333         {
1334             assert(!bInterMolInteractions);
1335             /* The vsite construction goes where the vsite itself is */
1336             if (iz == 0)
1337             {
1338                 add_vsite(dd->ga2la, index, rtil, ftype, nral,
1339                           TRUE, i, i_gl, i_mol,
1340                           iatoms, idef, vsite_pbc, vsite_pbc_nalloc);
1341             }
1342             j += 1 + nral + 2;
1343         }
1344         else
1345         {
1346             gmx_bool bUse;
1347
1348             /* Copy the type */
1349             tiatoms[0] = iatoms[0];
1350
1351             if (nral == 1)
1352             {
1353                 assert(!bInterMolInteractions);
1354                 /* Assign single-body interactions to the home zone */
1355                 if (iz == 0)
1356                 {
1357                     bUse       = TRUE;
1358                     tiatoms[1] = i;
1359                     if (ftype == F_POSRES)
1360                     {
1361                         add_posres(mol, i_mol, molb, tiatoms, ip_in,
1362                                    idef);
1363                     }
1364                     else if (ftype == F_FBPOSRES)
1365                     {
1366                         add_fbposres(mol, i_mol, molb, tiatoms, ip_in,
1367                                      idef);
1368                     }
1369                 }
1370                 else
1371                 {
1372                     bUse = FALSE;
1373                 }
1374             }
1375             else if (nral == 2)
1376             {
1377                 /* This is a two-body interaction, we can assign
1378                  * analogous to the non-bonded assignments.
1379                  */
1380                 int k_gl, a_loc, kz;
1381
1382                 if (!bInterMolInteractions)
1383                 {
1384                     /* Get the global index using the offset in the molecule */
1385                     k_gl = i_gl + iatoms[2] - i_mol;
1386                 }
1387                 else
1388                 {
1389                     k_gl = iatoms[2];
1390                 }
1391                 if (!ga2la_get(dd->ga2la, k_gl, &a_loc, &kz))
1392                 {
1393                     bUse = FALSE;
1394                 }
1395                 else
1396                 {
1397                     if (kz >= zones->n)
1398                     {
1399                         kz -= zones->n;
1400                     }
1401                     /* Check zone interaction assignments */
1402                     bUse = ((iz < zones->nizone &&
1403                              iz <= kz &&
1404                              kz >= zones->izone[iz].j0 &&
1405                              kz <  zones->izone[iz].j1) ||
1406                             (kz < zones->nizone &&
1407                                   iz > kz &&
1408                              iz >= zones->izone[kz].j0 &&
1409                              iz <  zones->izone[kz].j1));
1410                     if (bUse)
1411                     {
1412                         tiatoms[1] = i;
1413                         tiatoms[2] = a_loc;
1414                         /* If necessary check the cgcm distance */
1415                         if (bRCheck2B &&
1416                             dd_dist2(pbc_null, cg_cm, la2lc,
1417                                      tiatoms[1], tiatoms[2]) >= rc2)
1418                         {
1419                             bUse = FALSE;
1420                         }
1421                     }
1422                 }
1423             }
1424             else
1425             {
1426                 /* Assign this multi-body bonded interaction to
1427                  * the local node if we have all the atoms involved
1428                  * (local or communicated) and the minimum zone shift
1429                  * in each dimension is zero, for dimensions
1430                  * with 2 DD cells an extra check may be necessary.
1431                  */
1432                 ivec k_zero, k_plus;
1433                 int  k;
1434
1435                 bUse = TRUE;
1436                 clear_ivec(k_zero);
1437                 clear_ivec(k_plus);
1438                 for (k = 1; k <= nral && bUse; k++)
1439                 {
1440                     gmx_bool bLocal;
1441                     int      k_gl, a_loc;
1442                     int      kz;
1443
1444                     if (!bInterMolInteractions)
1445                     {
1446                         /* Get the global index using the offset in the molecule */
1447                         k_gl = i_gl + iatoms[k] - i_mol;
1448                     }
1449                     else
1450                     {
1451                         k_gl = iatoms[k];
1452                     }
1453                     bLocal = ga2la_get(dd->ga2la, k_gl, &a_loc, &kz);
1454                     if (!bLocal || kz >= zones->n)
1455                     {
1456                         /* We do not have this atom of this interaction
1457                          * locally, or it comes from more than one cell
1458                          * away.
1459                          */
1460                         bUse = FALSE;
1461                     }
1462                     else
1463                     {
1464                         int d;
1465
1466                         tiatoms[k] = a_loc;
1467                         for (d = 0; d < DIM; d++)
1468                         {
1469                             if (zones->shift[kz][d] == 0)
1470                             {
1471                                 k_zero[d] = k;
1472                             }
1473                             else
1474                             {
1475                                 k_plus[d] = k;
1476                             }
1477                         }
1478                     }
1479                 }
1480                 bUse = (bUse &&
1481                         k_zero[XX] && k_zero[YY] && k_zero[ZZ]);
1482                 if (bRCheckMB)
1483                 {
1484                     int d;
1485
1486                     for (d = 0; (d < DIM && bUse); d++)
1487                     {
1488                         /* Check if the cg_cm distance falls within
1489                          * the cut-off to avoid possible multiple
1490                          * assignments of bonded interactions.
1491                          */
1492                         if (rcheck[d] &&
1493                             k_plus[d] &&
1494                             dd_dist2(pbc_null, cg_cm, la2lc,
1495                                      tiatoms[k_zero[d]], tiatoms[k_plus[d]]) >= rc2)
1496                         {
1497                             bUse = FALSE;
1498                         }
1499                     }
1500                 }
1501             }
1502             if (bUse)
1503             {
1504                 /* Add this interaction to the local topology */
1505                 add_ifunc(nral, tiatoms, &idef->il[ftype]);
1506                 /* Sum so we can check in global_stat
1507                  * if we have everything.
1508                  */
1509                 if (bBCheck ||
1510                     !(interaction_function[ftype].flags & IF_LIMZERO))
1511                 {
1512                     (*nbonded_local)++;
1513                 }
1514             }
1515             j += 1 + nral;
1516         }
1517     }
1518 }
1519
1520 /*! \brief This function looks up and assigns bonded interactions for zone iz.
1521  *
1522  * With thread parallelizing each thread acts on a different atom range:
1523  * at_start to at_end.
1524  */
1525 static int make_bondeds_zone(gmx_domdec_t *dd,
1526                              const gmx_domdec_zones_t *zones,
1527                              const std::vector<gmx_molblock_t> &molb,
1528                              gmx_bool bRCheckMB, ivec rcheck, gmx_bool bRCheck2B,
1529                              real rc2,
1530                              int *la2lc, t_pbc *pbc_null, rvec *cg_cm,
1531                              const t_iparams *ip_in,
1532                              t_idef *idef,
1533                              int **vsite_pbc,
1534                              int *vsite_pbc_nalloc,
1535                              int izone,
1536                              int at_start, int at_end)
1537 {
1538     int                i, i_gl, mb, mt, mol, i_mol;
1539     int               *index, *rtil;
1540     gmx_bool           bBCheck;
1541     gmx_reverse_top_t *rt;
1542     int                nbonded_local;
1543
1544     rt = dd->reverse_top;
1545
1546     bBCheck = rt->bBCheck;
1547
1548     nbonded_local = 0;
1549
1550     for (i = at_start; i < at_end; i++)
1551     {
1552         /* Get the global atom number */
1553         i_gl = dd->gatindex[i];
1554         global_atomnr_to_moltype_ind(rt, i_gl, &mb, &mt, &mol, &i_mol);
1555         /* Check all intramolecular interactions assigned to this atom */
1556         index = rt->ril_mt[mt].index;
1557         rtil  = rt->ril_mt[mt].il;
1558
1559         check_assign_interactions_atom(i, i_gl, mol, i_mol,
1560                                        index, rtil, FALSE,
1561                                        index[i_mol], index[i_mol+1],
1562                                        dd, zones,
1563                                        &molb[mb],
1564                                        bRCheckMB, rcheck, bRCheck2B, rc2,
1565                                        la2lc,
1566                                        pbc_null,
1567                                        cg_cm,
1568                                        ip_in,
1569                                        idef, vsite_pbc, vsite_pbc_nalloc,
1570                                        izone,
1571                                        bBCheck,
1572                                        &nbonded_local);
1573
1574
1575         if (rt->bIntermolecularInteractions)
1576         {
1577             /* Check all intermolecular interactions assigned to this atom */
1578             index = rt->ril_intermol.index;
1579             rtil  = rt->ril_intermol.il;
1580
1581             check_assign_interactions_atom(i, i_gl, mol, i_mol,
1582                                            index, rtil, TRUE,
1583                                            index[i_gl], index[i_gl + 1],
1584                                            dd, zones,
1585                                            &molb[mb],
1586                                            bRCheckMB, rcheck, bRCheck2B, rc2,
1587                                            la2lc,
1588                                            pbc_null,
1589                                            cg_cm,
1590                                            ip_in,
1591                                            idef, vsite_pbc, vsite_pbc_nalloc,
1592                                            izone,
1593                                            bBCheck,
1594                                            &nbonded_local);
1595         }
1596     }
1597
1598     return nbonded_local;
1599 }
1600
1601 /*! \brief Set the exclusion data for i-zone \p iz for the case of no exclusions */
1602 static void set_no_exclusions_zone(gmx_domdec_t *dd, gmx_domdec_zones_t *zones,
1603                                    int iz, t_blocka *lexcls)
1604 {
1605     int  a0, a1, a;
1606
1607     a0 = dd->cgindex[zones->cg_range[iz]];
1608     a1 = dd->cgindex[zones->cg_range[iz+1]];
1609
1610     for (a = a0+1; a < a1+1; a++)
1611     {
1612         lexcls->index[a] = lexcls->nra;
1613     }
1614 }
1615
1616 /*! \brief Set the exclusion data for i-zone \p iz
1617  *
1618  * This is a legacy version for the group scheme of the same routine below.
1619  * Here charge groups and distance checks to ensure unique exclusions
1620  * are supported.
1621  */
1622 static int make_exclusions_zone_cg(gmx_domdec_t *dd, gmx_domdec_zones_t *zones,
1623                                    const std::vector<gmx_moltype_t> &moltype,
1624                                    gmx_bool bRCheck, real rc2,
1625                                    int *la2lc, t_pbc *pbc_null, rvec *cg_cm,
1626                                    const int *cginfo,
1627                                    t_blocka *lexcls,
1628                                    int iz,
1629                                    int cg_start, int cg_end)
1630 {
1631     int               n_excl_at_max, n, count, jla0, jla1, jla;
1632     int               cg, la0, la1, la, a_gl, mb, mt, mol, a_mol, j, aj_mol;
1633     const t_blocka   *excls;
1634     gmx_ga2la_t      *ga2la;
1635     int               cell;
1636
1637     ga2la = dd->ga2la;
1638
1639     jla0 = dd->cgindex[zones->izone[iz].jcg0];
1640     jla1 = dd->cgindex[zones->izone[iz].jcg1];
1641
1642     n_excl_at_max = dd->reverse_top->n_excl_at_max;
1643
1644     /* We set the end index, but note that we might not start at zero here */
1645     lexcls->nr = dd->cgindex[cg_end];
1646
1647     n     = lexcls->nra;
1648     count = 0;
1649     for (cg = cg_start; cg < cg_end; cg++)
1650     {
1651         if (n + (cg_end - cg_start)*n_excl_at_max > lexcls->nalloc_a)
1652         {
1653             lexcls->nalloc_a = over_alloc_large(n + (cg_end - cg_start)*n_excl_at_max);
1654             srenew(lexcls->a, lexcls->nalloc_a);
1655         }
1656         la0 = dd->cgindex[cg];
1657         la1 = dd->cgindex[cg+1];
1658         if (GET_CGINFO_EXCL_INTER(cginfo[cg]) ||
1659             !GET_CGINFO_EXCL_INTRA(cginfo[cg]))
1660         {
1661             /* Copy the exclusions from the global top */
1662             for (la = la0; la < la1; la++)
1663             {
1664                 lexcls->index[la] = n;
1665                 a_gl              = dd->gatindex[la];
1666                 global_atomnr_to_moltype_ind(dd->reverse_top, a_gl, &mb, &mt, &mol, &a_mol);
1667                 excls = &moltype[mt].excls;
1668                 for (j = excls->index[a_mol]; j < excls->index[a_mol+1]; j++)
1669                 {
1670                     aj_mol = excls->a[j];
1671                     /* This computation of jla is only correct intra-cg */
1672                     jla = la + aj_mol - a_mol;
1673                     if (jla >= la0 && jla < la1)
1674                     {
1675                         /* This is an intra-cg exclusion. We can skip
1676                          *  the global indexing and distance checking.
1677                          */
1678                         /* Intra-cg exclusions are only required
1679                          * for the home zone.
1680                          */
1681                         if (iz == 0)
1682                         {
1683                             lexcls->a[n++] = jla;
1684                             /* Check to avoid double counts */
1685                             if (jla > la)
1686                             {
1687                                 count++;
1688                             }
1689                         }
1690                     }
1691                     else
1692                     {
1693                         /* This is a inter-cg exclusion */
1694                         /* Since exclusions are pair interactions,
1695                          * just like non-bonded interactions,
1696                          * they can be assigned properly up
1697                          * to the DD cutoff (not cutoff_min as
1698                          * for the other bonded interactions).
1699                          */
1700                         if (ga2la_get(ga2la, a_gl+aj_mol-a_mol, &jla, &cell))
1701                         {
1702                             if (iz == 0 && cell == 0)
1703                             {
1704                                 lexcls->a[n++] = jla;
1705                                 /* Check to avoid double counts */
1706                                 if (jla > la)
1707                                 {
1708                                     count++;
1709                                 }
1710                             }
1711                             else if (jla >= jla0 && jla < jla1 &&
1712                                      (!bRCheck ||
1713                                       dd_dist2(pbc_null, cg_cm, la2lc, la, jla) < rc2))
1714                             {
1715                                 /* jla > la, since jla0 > la */
1716                                 lexcls->a[n++] = jla;
1717                                 count++;
1718                             }
1719                         }
1720                     }
1721                 }
1722             }
1723         }
1724         else
1725         {
1726             /* There are no inter-cg excls and this cg is self-excluded.
1727              * These exclusions are only required for zone 0,
1728              * since other zones do not see themselves.
1729              */
1730             if (iz == 0)
1731             {
1732                 for (la = la0; la < la1; la++)
1733                 {
1734                     lexcls->index[la] = n;
1735                     for (j = la0; j < la1; j++)
1736                     {
1737                         lexcls->a[n++] = j;
1738                     }
1739                 }
1740                 count += ((la1 - la0)*(la1 - la0 - 1))/2;
1741             }
1742             else
1743             {
1744                 /* We don't need exclusions for this cg */
1745                 for (la = la0; la < la1; la++)
1746                 {
1747                     lexcls->index[la] = n;
1748                 }
1749             }
1750         }
1751     }
1752
1753     lexcls->index[lexcls->nr] = n;
1754     lexcls->nra               = n;
1755
1756     return count;
1757 }
1758
1759 /*! \brief Set the exclusion data for i-zone \p iz */
1760 static void make_exclusions_zone(gmx_domdec_t *dd,
1761                                  gmx_domdec_zones_t *zones,
1762                                  const std::vector<gmx_moltype_t> &moltype,
1763                                  const int *cginfo,
1764                                  t_blocka *lexcls,
1765                                  int iz,
1766                                  int at_start, int at_end)
1767 {
1768     gmx_ga2la_t *ga2la;
1769     int          jla0, jla1;
1770     int          n_excl_at_max, n, at;
1771
1772     ga2la = dd->ga2la;
1773
1774     jla0 = dd->cgindex[zones->izone[iz].jcg0];
1775     jla1 = dd->cgindex[zones->izone[iz].jcg1];
1776
1777     n_excl_at_max = dd->reverse_top->n_excl_at_max;
1778
1779     /* We set the end index, but note that we might not start at zero here */
1780     lexcls->nr = at_end;
1781
1782     n = lexcls->nra;
1783     for (at = at_start; at < at_end; at++)
1784     {
1785         if (n + 1000 > lexcls->nalloc_a)
1786         {
1787             lexcls->nalloc_a = over_alloc_large(n + 1000);
1788             srenew(lexcls->a, lexcls->nalloc_a);
1789         }
1790         if (GET_CGINFO_EXCL_INTER(cginfo[at]))
1791         {
1792             int             a_gl, mb, mt, mol, a_mol, j;
1793             const t_blocka *excls;
1794
1795             if (n + n_excl_at_max > lexcls->nalloc_a)
1796             {
1797                 lexcls->nalloc_a = over_alloc_large(n + n_excl_at_max);
1798                 srenew(lexcls->a, lexcls->nalloc_a);
1799             }
1800
1801             /* Copy the exclusions from the global top */
1802             lexcls->index[at] = n;
1803             a_gl              = dd->gatindex[at];
1804             global_atomnr_to_moltype_ind(dd->reverse_top, a_gl,
1805                                          &mb, &mt, &mol, &a_mol);
1806             excls = &moltype[mt].excls;
1807             for (j = excls->index[a_mol]; j < excls->index[a_mol + 1]; j++)
1808             {
1809                 int aj_mol, at_j, cell;
1810
1811                 aj_mol = excls->a[j];
1812
1813                 if (ga2la_get(ga2la, a_gl + aj_mol - a_mol, &at_j, &cell))
1814                 {
1815                     /* This check is not necessary, but it can reduce
1816                      * the number of exclusions in the list, which in turn
1817                      * can speed up the pair list construction a bit.
1818                      */
1819                     if (at_j >= jla0 && at_j < jla1)
1820                     {
1821                         lexcls->a[n++] = at_j;
1822                     }
1823                 }
1824             }
1825         }
1826         else
1827         {
1828             /* We don't need exclusions for this atom */
1829             lexcls->index[at] = n;
1830         }
1831     }
1832
1833     lexcls->index[lexcls->nr] = n;
1834     lexcls->nra               = n;
1835 }
1836
1837
1838 /*! \brief Ensure we have enough space in \p ba for \p nindex_max indices */
1839 static void check_alloc_index(t_blocka *ba, int nindex_max)
1840 {
1841     if (nindex_max+1 > ba->nalloc_index)
1842     {
1843         ba->nalloc_index = over_alloc_dd(nindex_max+1);
1844         srenew(ba->index, ba->nalloc_index);
1845     }
1846 }
1847
1848 /*! \brief Ensure that we have enough space for exclusion storate in \p lexcls */
1849 static void check_exclusions_alloc(gmx_domdec_t *dd, gmx_domdec_zones_t *zones,
1850                                    t_blocka *lexcls)
1851 {
1852     int nr;
1853     int thread;
1854
1855     nr = dd->cgindex[zones->izone[zones->nizone-1].cg1];
1856
1857     check_alloc_index(lexcls, nr);
1858
1859     for (thread = 1; thread < dd->reverse_top->nthread; thread++)
1860     {
1861         check_alloc_index(&dd->reverse_top->th_work[thread].excl, nr);
1862     }
1863 }
1864
1865 /*! \brief Set the total count indexes for the local exclusions, needed by several functions */
1866 static void finish_local_exclusions(gmx_domdec_t *dd, gmx_domdec_zones_t *zones,
1867                                     t_blocka *lexcls)
1868 {
1869     int la0, la;
1870
1871     lexcls->nr = dd->cgindex[zones->izone[zones->nizone-1].cg1];
1872
1873     if (dd->n_intercg_excl == 0)
1874     {
1875         /* There are no exclusions involving non-home charge groups,
1876          * but we need to set the indices for neighborsearching.
1877          */
1878         la0 = dd->cgindex[zones->izone[0].cg1];
1879         for (la = la0; la < lexcls->nr; la++)
1880         {
1881             lexcls->index[la] = lexcls->nra;
1882         }
1883
1884         /* nr is only used to loop over the exclusions for Ewald and RF,
1885          * so we can set it to the number of home atoms for efficiency.
1886          */
1887         lexcls->nr = dd->cgindex[zones->izone[0].cg1];
1888     }
1889 }
1890
1891 /*! \brief Clear a t_idef struct */
1892 static void clear_idef(t_idef *idef)
1893 {
1894     int  ftype;
1895
1896     /* Clear the counts */
1897     for (ftype = 0; ftype < F_NRE; ftype++)
1898     {
1899         idef->il[ftype].nr = 0;
1900     }
1901 }
1902
1903 /*! \brief Generate and store all required local bonded interactions in \p idef and local exclusions in \p lexcls */
1904 static int make_local_bondeds_excls(gmx_domdec_t *dd,
1905                                     gmx_domdec_zones_t *zones,
1906                                     const gmx_mtop_t *mtop,
1907                                     const int *cginfo,
1908                                     gmx_bool bRCheckMB, ivec rcheck, gmx_bool bRCheck2B,
1909                                     real rc,
1910                                     int *la2lc, t_pbc *pbc_null, rvec *cg_cm,
1911                                     t_idef *idef, gmx_vsite_t *vsite,
1912                                     t_blocka *lexcls, int *excl_count)
1913 {
1914     int                nzone_bondeds, nzone_excl;
1915     int                izone, cg0, cg1;
1916     real               rc2;
1917     int                nbonded_local;
1918     int                thread;
1919     gmx_reverse_top_t *rt;
1920
1921     if (dd->reverse_top->bInterCGInteractions)
1922     {
1923         nzone_bondeds = zones->n;
1924     }
1925     else
1926     {
1927         /* Only single charge group (or atom) molecules, so interactions don't
1928          * cross zone boundaries and we only need to assign in the home zone.
1929          */
1930         nzone_bondeds = 1;
1931     }
1932
1933     if (dd->n_intercg_excl > 0)
1934     {
1935         /* We only use exclusions from i-zones to i- and j-zones */
1936         nzone_excl = zones->nizone;
1937     }
1938     else
1939     {
1940         /* There are no inter-cg exclusions and only zone 0 sees itself */
1941         nzone_excl = 1;
1942     }
1943
1944     check_exclusions_alloc(dd, zones, lexcls);
1945
1946     rt = dd->reverse_top;
1947
1948     rc2 = rc*rc;
1949
1950     /* Clear the counts */
1951     clear_idef(idef);
1952     nbonded_local = 0;
1953
1954     lexcls->nr    = 0;
1955     lexcls->nra   = 0;
1956     *excl_count   = 0;
1957
1958     for (izone = 0; izone < nzone_bondeds; izone++)
1959     {
1960         cg0 = zones->cg_range[izone];
1961         cg1 = zones->cg_range[izone + 1];
1962
1963 #pragma omp parallel for num_threads(rt->nthread) schedule(static)
1964         for (thread = 0; thread < rt->nthread; thread++)
1965         {
1966             try
1967             {
1968                 int       cg0t, cg1t;
1969                 t_idef   *idef_t;
1970                 int     **vsite_pbc;
1971                 int      *vsite_pbc_nalloc;
1972                 t_blocka *excl_t;
1973
1974                 cg0t = cg0 + ((cg1 - cg0)* thread   )/rt->nthread;
1975                 cg1t = cg0 + ((cg1 - cg0)*(thread+1))/rt->nthread;
1976
1977                 if (thread == 0)
1978                 {
1979                     idef_t = idef;
1980                 }
1981                 else
1982                 {
1983                     idef_t = &rt->th_work[thread].idef;
1984                     clear_idef(idef_t);
1985                 }
1986
1987                 if (vsite && vsite->bHaveChargeGroups && vsite->n_intercg_vsite > 0)
1988                 {
1989                     if (thread == 0)
1990                     {
1991                         vsite_pbc        = vsite->vsite_pbc_loc;
1992                         vsite_pbc_nalloc = vsite->vsite_pbc_loc_nalloc;
1993                     }
1994                     else
1995                     {
1996                         vsite_pbc        = rt->th_work[thread].vsite_pbc;
1997                         vsite_pbc_nalloc = rt->th_work[thread].vsite_pbc_nalloc;
1998                     }
1999                 }
2000                 else
2001                 {
2002                     vsite_pbc        = nullptr;
2003                     vsite_pbc_nalloc = nullptr;
2004                 }
2005
2006                 rt->th_work[thread].nbonded =
2007                     make_bondeds_zone(dd, zones,
2008                                       mtop->molblock,
2009                                       bRCheckMB, rcheck, bRCheck2B, rc2,
2010                                       la2lc, pbc_null, cg_cm, idef->iparams,
2011                                       idef_t,
2012                                       vsite_pbc, vsite_pbc_nalloc,
2013                                       izone,
2014                                       dd->cgindex[cg0t], dd->cgindex[cg1t]);
2015
2016                 if (izone < nzone_excl)
2017                 {
2018                     if (thread == 0)
2019                     {
2020                         excl_t = lexcls;
2021                     }
2022                     else
2023                     {
2024                         excl_t      = &rt->th_work[thread].excl;
2025                         excl_t->nr  = 0;
2026                         excl_t->nra = 0;
2027                     }
2028
2029                     if (dd->cgindex[dd->ncg_tot] == dd->ncg_tot &&
2030                         !rt->bExclRequired)
2031                     {
2032                         /* No charge groups and no distance check required */
2033                         make_exclusions_zone(dd, zones,
2034                                              mtop->moltype, cginfo,
2035                                              excl_t,
2036                                              izone,
2037                                              cg0t, cg1t);
2038                     }
2039                     else
2040                     {
2041                         rt->th_work[thread].excl_count =
2042                             make_exclusions_zone_cg(dd, zones,
2043                                                     mtop->moltype, bRCheck2B, rc2,
2044                                                     la2lc, pbc_null, cg_cm, cginfo,
2045                                                     excl_t,
2046                                                     izone,
2047                                                     cg0t, cg1t);
2048                     }
2049                 }
2050             }
2051             GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
2052         }
2053
2054         if (rt->nthread > 1)
2055         {
2056             combine_idef(idef, rt->th_work, rt->nthread, vsite);
2057         }
2058
2059         for (thread = 0; thread < rt->nthread; thread++)
2060         {
2061             nbonded_local += rt->th_work[thread].nbonded;
2062         }
2063
2064         if (izone < nzone_excl)
2065         {
2066             if (rt->nthread > 1)
2067             {
2068                 combine_blocka(lexcls, rt->th_work, rt->nthread);
2069             }
2070
2071             for (thread = 0; thread < rt->nthread; thread++)
2072             {
2073                 *excl_count += rt->th_work[thread].excl_count;
2074             }
2075         }
2076     }
2077
2078     /* Some zones might not have exclusions, but some code still needs to
2079      * loop over the index, so we set the indices here.
2080      */
2081     for (izone = nzone_excl; izone < zones->nizone; izone++)
2082     {
2083         set_no_exclusions_zone(dd, zones, izone, lexcls);
2084     }
2085
2086     finish_local_exclusions(dd, zones, lexcls);
2087     if (debug)
2088     {
2089         fprintf(debug, "We have %d exclusions, check count %d\n",
2090                 lexcls->nra, *excl_count);
2091     }
2092
2093     return nbonded_local;
2094 }
2095
2096 void dd_make_local_cgs(gmx_domdec_t *dd, t_block *lcgs)
2097 {
2098     lcgs->nr    = dd->ncg_tot;
2099     lcgs->index = dd->cgindex;
2100 }
2101
2102 void dd_make_local_top(gmx_domdec_t *dd, gmx_domdec_zones_t *zones,
2103                        int npbcdim, matrix box,
2104                        rvec cellsize_min, ivec npulse,
2105                        t_forcerec *fr,
2106                        rvec *cgcm_or_x,
2107                        gmx_vsite_t *vsite,
2108                        const gmx_mtop_t *mtop, gmx_localtop_t *ltop)
2109 {
2110     gmx_bool bRCheckMB, bRCheck2B;
2111     real     rc = -1;
2112     ivec     rcheck;
2113     int      d, nexcl;
2114     t_pbc    pbc, *pbc_null = nullptr;
2115
2116     if (debug)
2117     {
2118         fprintf(debug, "Making local topology\n");
2119     }
2120
2121     dd_make_local_cgs(dd, &ltop->cgs);
2122
2123     bRCheckMB   = FALSE;
2124     bRCheck2B   = FALSE;
2125
2126     if (dd->reverse_top->bInterCGInteractions)
2127     {
2128         /* We need to check to which cell bondeds should be assigned */
2129         rc = dd_cutoff_twobody(dd);
2130         if (debug)
2131         {
2132             fprintf(debug, "Two-body bonded cut-off distance is %g\n", rc);
2133         }
2134
2135         /* Should we check cg_cm distances when assigning bonded interactions? */
2136         for (d = 0; d < DIM; d++)
2137         {
2138             rcheck[d] = FALSE;
2139             /* Only need to check for dimensions where the part of the box
2140              * that is not communicated is smaller than the cut-off.
2141              */
2142             if (d < npbcdim && dd->nc[d] > 1 &&
2143                 (dd->nc[d] - npulse[d])*cellsize_min[d] < 2*rc)
2144             {
2145                 if (dd->nc[d] == 2)
2146                 {
2147                     rcheck[d] = TRUE;
2148                     bRCheckMB = TRUE;
2149                 }
2150                 /* Check for interactions between two atoms,
2151                  * where we can allow interactions up to the cut-off,
2152                  * instead of up to the smallest cell dimension.
2153                  */
2154                 bRCheck2B = TRUE;
2155             }
2156             if (debug)
2157             {
2158                 fprintf(debug,
2159                         "dim %d cellmin %f bonded rcheck[%d] = %d, bRCheck2B = %d\n",
2160                         d, cellsize_min[d], d, rcheck[d], bRCheck2B);
2161             }
2162         }
2163         if (bRCheckMB || bRCheck2B)
2164         {
2165             make_la2lc(dd);
2166             if (fr->bMolPBC)
2167             {
2168                 pbc_null = set_pbc_dd(&pbc, fr->ePBC, dd->nc, TRUE, box);
2169             }
2170             else
2171             {
2172                 pbc_null = nullptr;
2173             }
2174         }
2175     }
2176
2177     dd->nbonded_local =
2178         make_local_bondeds_excls(dd, zones, mtop, fr->cginfo,
2179                                  bRCheckMB, rcheck, bRCheck2B, rc,
2180                                  dd->la2lc,
2181                                  pbc_null, cgcm_or_x,
2182                                  &ltop->idef, vsite,
2183                                  &ltop->excls, &nexcl);
2184
2185     /* The ilist is not sorted yet,
2186      * we can only do this when we have the charge arrays.
2187      */
2188     ltop->idef.ilsort = ilsortUNKNOWN;
2189
2190     if (dd->reverse_top->bExclRequired)
2191     {
2192         dd->nbonded_local += nexcl;
2193     }
2194
2195     ltop->atomtypes  = mtop->atomtypes;
2196 }
2197
2198 void dd_sort_local_top(gmx_domdec_t *dd, const t_mdatoms *mdatoms,
2199                        gmx_localtop_t *ltop)
2200 {
2201     if (dd->reverse_top->ilsort == ilsortNO_FE)
2202     {
2203         ltop->idef.ilsort = ilsortNO_FE;
2204     }
2205     else
2206     {
2207         gmx_sort_ilist_fe(&ltop->idef, mdatoms->chargeA, mdatoms->chargeB);
2208     }
2209 }
2210
2211 gmx_localtop_t *dd_init_local_top(const gmx_mtop_t *top_global)
2212 {
2213     gmx_localtop_t *top;
2214     int             i;
2215
2216     snew(top, 1);
2217
2218     top->idef.ntypes    = top_global->ffparams.ntypes;
2219     top->idef.atnr      = top_global->ffparams.atnr;
2220     top->idef.functype  = top_global->ffparams.functype;
2221     top->idef.iparams   = top_global->ffparams.iparams;
2222     top->idef.fudgeQQ   = top_global->ffparams.fudgeQQ;
2223     top->idef.cmap_grid = top_global->ffparams.cmap_grid;
2224
2225     for (i = 0; i < F_NRE; i++)
2226     {
2227         top->idef.il[i].iatoms = nullptr;
2228         top->idef.il[i].nalloc = 0;
2229     }
2230     top->idef.ilsort   = ilsortUNKNOWN;
2231
2232     return top;
2233 }
2234
2235 void dd_init_local_state(gmx_domdec_t *dd,
2236                          t_state *state_global, t_state *state_local)
2237 {
2238     int buf[NITEM_DD_INIT_LOCAL_STATE];
2239
2240     if (DDMASTER(dd))
2241     {
2242         buf[0] = state_global->flags;
2243         buf[1] = state_global->ngtc;
2244         buf[2] = state_global->nnhpres;
2245         buf[3] = state_global->nhchainlength;
2246         buf[4] = state_global->dfhist ? state_global->dfhist->nlambda : 0;
2247     }
2248     dd_bcast(dd, NITEM_DD_INIT_LOCAL_STATE*sizeof(int), buf);
2249
2250     init_gtc_state(state_local, buf[1], buf[2], buf[3]);
2251     init_dfhist_state(state_local, buf[4]);
2252     state_local->flags = buf[0];
2253 }
2254
2255 /*! \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 */
2256 static void check_link(t_blocka *link, int cg_gl, int cg_gl_j)
2257 {
2258     int      k;
2259     gmx_bool bFound;
2260
2261     bFound = FALSE;
2262     for (k = link->index[cg_gl]; k < link->index[cg_gl+1]; k++)
2263     {
2264         GMX_RELEASE_ASSERT(link->a, "Inconsistent NULL pointer while making charge-group links");
2265         if (link->a[k] == cg_gl_j)
2266         {
2267             bFound = TRUE;
2268         }
2269     }
2270     if (!bFound)
2271     {
2272         GMX_RELEASE_ASSERT(link->a || link->index[cg_gl+1]+1 > link->nalloc_a,
2273                            "Inconsistent allocation of link");
2274         /* Add this charge group link */
2275         if (link->index[cg_gl+1]+1 > link->nalloc_a)
2276         {
2277             link->nalloc_a = over_alloc_large(link->index[cg_gl+1]+1);
2278             srenew(link->a, link->nalloc_a);
2279         }
2280         link->a[link->index[cg_gl+1]] = cg_gl_j;
2281         link->index[cg_gl+1]++;
2282     }
2283 }
2284
2285 /*! \brief Return a vector of the charge group index for all atoms */
2286 static std::vector<int> make_at2cg(const t_block &cgs)
2287 {
2288     std::vector<int> at2cg(cgs.index[cgs.nr]);
2289     for (int cg = 0; cg < cgs.nr; cg++)
2290     {
2291         for (int a = cgs.index[cg]; a < cgs.index[cg + 1]; a++)
2292         {
2293             at2cg[a] = cg;
2294         }
2295     }
2296
2297     return at2cg;
2298 }
2299
2300 t_blocka *make_charge_group_links(const gmx_mtop_t *mtop, gmx_domdec_t *dd,
2301                                   cginfo_mb_t *cginfo_mb)
2302 {
2303     gmx_bool            bExclRequired;
2304     reverse_ilist_t     ril, ril_intermol;
2305     t_blocka           *link;
2306     cginfo_mb_t        *cgi_mb;
2307
2308     /* For each charge group make a list of other charge groups
2309      * in the system that a linked to it via bonded interactions
2310      * which are also stored in reverse_top.
2311      */
2312
2313     bExclRequired = dd->reverse_top->bExclRequired;
2314
2315     if (mtop->bIntermolecularInteractions)
2316     {
2317         if (ncg_mtop(mtop) < mtop->natoms)
2318         {
2319             gmx_fatal(FARGS, "The combination of intermolecular interactions, charge groups and domain decomposition is not supported. Use cutoff-scheme=Verlet (which removes the charge groups) or run without domain decomposition.");
2320         }
2321
2322         t_atoms atoms;
2323
2324         atoms.nr   = mtop->natoms;
2325         atoms.atom = nullptr;
2326
2327         make_reverse_ilist(mtop->intermolecular_ilist, &atoms,
2328                            nullptr, FALSE, FALSE, FALSE, TRUE, &ril_intermol);
2329     }
2330
2331     snew(link, 1);
2332     snew(link->index, ncg_mtop(mtop)+1);
2333     link->nalloc_a = 0;
2334     link->a        = nullptr;
2335
2336     link->index[0] = 0;
2337     int cg_offset  = 0;
2338     int ncgi       = 0;
2339     for (size_t mb = 0; mb < mtop->molblock.size(); mb++)
2340     {
2341         const gmx_molblock_t &molb = mtop->molblock[mb];
2342         if (molb.nmol == 0)
2343         {
2344             continue;
2345         }
2346         const gmx_moltype_t &molt  = mtop->moltype[molb.type];
2347         const t_block       &cgs   = molt.cgs;
2348         const t_blocka      &excls = molt.excls;
2349         std::vector<int>     a2c   = make_at2cg(cgs);
2350         /* Make a reverse ilist in which the interactions are linked
2351          * to all atoms, not only the first atom as in gmx_reverse_top.
2352          * The constraints are discarded here.
2353          */
2354         make_reverse_ilist(molt.ilist, &molt.atoms,
2355                            nullptr, FALSE, FALSE, FALSE, TRUE, &ril);
2356
2357         cgi_mb = &cginfo_mb[mb];
2358
2359         int mol;
2360         for (mol = 0; mol < (mtop->bIntermolecularInteractions ? molb.nmol : 1); mol++)
2361         {
2362             for (int cg = 0; cg < cgs.nr; cg++)
2363             {
2364                 int cg_gl            = cg_offset + cg;
2365                 link->index[cg_gl+1] = link->index[cg_gl];
2366                 for (int a = cgs.index[cg]; a < cgs.index[cg + 1]; a++)
2367                 {
2368                     int i = ril.index[a];
2369                     while (i < ril.index[a+1])
2370                     {
2371                         int ftype = ril.il[i++];
2372                         int nral  = NRAL(ftype);
2373                         /* Skip the ifunc index */
2374                         i++;
2375                         for (int j = 0; j < nral; j++)
2376                         {
2377                             int aj = ril.il[i + j];
2378                             if (a2c[aj] != cg)
2379                             {
2380                                 check_link(link, cg_gl, cg_offset+a2c[aj]);
2381                             }
2382                         }
2383                         i += nral_rt(ftype);
2384                     }
2385                     if (bExclRequired)
2386                     {
2387                         /* Exclusions always go both ways */
2388                         for (int j = excls.index[a]; j < excls.index[a + 1]; j++)
2389                         {
2390                             int aj = excls.a[j];
2391                             if (a2c[aj] != cg)
2392                             {
2393                                 check_link(link, cg_gl, cg_offset+a2c[aj]);
2394                             }
2395                         }
2396                     }
2397
2398                     if (mtop->bIntermolecularInteractions)
2399                     {
2400                         int i = ril_intermol.index[a];
2401                         while (i < ril_intermol.index[a+1])
2402                         {
2403                             int ftype = ril_intermol.il[i++];
2404                             int nral  = NRAL(ftype);
2405                             /* Skip the ifunc index */
2406                             i++;
2407                             for (int j = 0; j < nral; j++)
2408                             {
2409                                 /* Here we assume we have no charge groups;
2410                                  * this has been checked above.
2411                                  */
2412                                 int aj = ril_intermol.il[i + j];
2413                                 check_link(link, cg_gl, aj);
2414                             }
2415                             i += nral_rt(ftype);
2416                         }
2417                     }
2418                 }
2419                 if (link->index[cg_gl+1] - link->index[cg_gl] > 0)
2420                 {
2421                     SET_CGINFO_BOND_INTER(cgi_mb->cginfo[cg]);
2422                     ncgi++;
2423                 }
2424             }
2425
2426             cg_offset += cgs.nr;
2427         }
2428         int nlink_mol = link->index[cg_offset] - link->index[cg_offset - cgs.nr];
2429
2430         destroy_reverse_ilist(&ril);
2431
2432         if (debug)
2433         {
2434             fprintf(debug, "molecule type '%s' %d cgs has %d cg links through bonded interac.\n", *molt.name, cgs.nr, nlink_mol);
2435         }
2436
2437         if (molb.nmol > mol)
2438         {
2439             /* Copy the data for the rest of the molecules in this block */
2440             link->nalloc_a += (molb.nmol - mol)*nlink_mol;
2441             srenew(link->a, link->nalloc_a);
2442             for (; mol < molb.nmol; mol++)
2443             {
2444                 for (int cg = 0; cg < cgs.nr; cg++)
2445                 {
2446                     int cg_gl              = cg_offset + cg;
2447                     link->index[cg_gl + 1] =
2448                         link->index[cg_gl + 1 - cgs.nr] + nlink_mol;
2449                     for (int j = link->index[cg_gl]; j < link->index[cg_gl+1]; j++)
2450                     {
2451                         link->a[j] = link->a[j - nlink_mol] + cgs.nr;
2452                     }
2453                     if (link->index[cg_gl+1] - link->index[cg_gl] > 0 &&
2454                         cg_gl - cgi_mb->cg_start < cgi_mb->cg_mod)
2455                     {
2456                         SET_CGINFO_BOND_INTER(cgi_mb->cginfo[cg_gl - cgi_mb->cg_start]);
2457                         ncgi++;
2458                     }
2459                 }
2460                 cg_offset += cgs.nr;
2461             }
2462         }
2463     }
2464
2465     if (mtop->bIntermolecularInteractions)
2466     {
2467         destroy_reverse_ilist(&ril_intermol);
2468     }
2469
2470     if (debug)
2471     {
2472         fprintf(debug, "Of the %d charge groups %d are linked via bonded interactions\n", ncg_mtop(mtop), ncgi);
2473     }
2474
2475     return link;
2476 }
2477
2478 typedef struct {
2479     real r2;
2480     int  ftype;
2481     int  a1;
2482     int  a2;
2483 } bonded_distance_t;
2484
2485 /*! \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 */
2486 static void update_max_bonded_distance(real r2, int ftype, int a1, int a2,
2487                                        bonded_distance_t *bd)
2488 {
2489     if (r2 > bd->r2)
2490     {
2491         bd->r2    = r2;
2492         bd->ftype = ftype;
2493         bd->a1    = a1;
2494         bd->a2    = a2;
2495     }
2496 }
2497
2498 /*! \brief Set the distance, function type and atom indices for the longest distance between charge-groups of molecule type \p molt for two-body and multi-body bonded interactions */
2499 static void bonded_cg_distance_mol(const gmx_moltype_t *molt,
2500                                    const std::vector<int> &at2cg,
2501                                    gmx_bool bBCheck, gmx_bool bExcl, rvec *cg_cm,
2502                                    bonded_distance_t *bd_2b,
2503                                    bonded_distance_t *bd_mb)
2504 {
2505     for (int ftype = 0; ftype < F_NRE; ftype++)
2506     {
2507         if (dd_check_ftype(ftype, bBCheck, FALSE, FALSE))
2508         {
2509             const t_ilist *il   = &molt->ilist[ftype];
2510             int            nral = NRAL(ftype);
2511             if (nral > 1)
2512             {
2513                 for (int i = 0; i < il->nr; i += 1+nral)
2514                 {
2515                     for (int ai = 0; ai < nral; ai++)
2516                     {
2517                         int cgi = at2cg[il->iatoms[i+1+ai]];
2518                         for (int aj = ai + 1; aj < nral; aj++)
2519                         {
2520                             int cgj = at2cg[il->iatoms[i+1+aj]];
2521                             if (cgi != cgj)
2522                             {
2523                                 real rij2 = distance2(cg_cm[cgi], cg_cm[cgj]);
2524
2525                                 update_max_bonded_distance(rij2, ftype,
2526                                                            il->iatoms[i+1+ai],
2527                                                            il->iatoms[i+1+aj],
2528                                                            (nral == 2) ? bd_2b : bd_mb);
2529                             }
2530                         }
2531                     }
2532                 }
2533             }
2534         }
2535     }
2536     if (bExcl)
2537     {
2538         const t_blocka *excls = &molt->excls;
2539         for (int ai = 0; ai < excls->nr; ai++)
2540         {
2541             int cgi = at2cg[ai];
2542             for (int j = excls->index[ai]; j < excls->index[ai+1]; j++)
2543             {
2544                 int cgj = at2cg[excls->a[j]];
2545                 if (cgi != cgj)
2546                 {
2547                     real rij2 = distance2(cg_cm[cgi], cg_cm[cgj]);
2548
2549                     /* There is no function type for exclusions, use -1 */
2550                     update_max_bonded_distance(rij2, -1, ai, excls->a[j], bd_2b);
2551                 }
2552             }
2553         }
2554     }
2555 }
2556
2557 /*! \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 */
2558 static void bonded_distance_intermol(const t_ilist *ilists_intermol,
2559                                      gmx_bool bBCheck,
2560                                      const rvec *x, int ePBC, const matrix box,
2561                                      bonded_distance_t *bd_2b,
2562                                      bonded_distance_t *bd_mb)
2563 {
2564     t_pbc pbc;
2565
2566     set_pbc(&pbc, ePBC, box);
2567
2568     for (int ftype = 0; ftype < F_NRE; ftype++)
2569     {
2570         if (dd_check_ftype(ftype, bBCheck, FALSE, FALSE))
2571         {
2572             const t_ilist *il   = &ilists_intermol[ftype];
2573             int            nral = NRAL(ftype);
2574
2575             /* No nral>1 check here, since intermol interactions always
2576              * have nral>=2 (and the code is also correct for nral=1).
2577              */
2578             for (int i = 0; i < il->nr; i += 1+nral)
2579             {
2580                 for (int ai = 0; ai < nral; ai++)
2581                 {
2582                     int atom_i = il->iatoms[i + 1 + ai];
2583
2584                     for (int aj = ai + 1; aj < nral; aj++)
2585                     {
2586                         rvec dx;
2587                         real rij2;
2588
2589                         int  atom_j = il->iatoms[i + 1 + aj];
2590
2591                         pbc_dx(&pbc, x[atom_i], x[atom_j], dx);
2592
2593                         rij2 = norm2(dx);
2594
2595                         update_max_bonded_distance(rij2, ftype,
2596                                                    atom_i, atom_j,
2597                                                    (nral == 2) ? bd_2b : bd_mb);
2598                     }
2599                 }
2600             }
2601         }
2602     }
2603 }
2604
2605 //! Returns whether \p molt has at least one virtual site
2606 static bool moltypeHasVsite(const gmx_moltype_t &molt)
2607 {
2608     bool hasVsite = false;
2609     for (int i = 0; i < F_NRE; i++)
2610     {
2611         if ((interaction_function[i].flags & IF_VSITE) &&
2612             molt.ilist[i].nr > 0)
2613         {
2614             hasVsite = true;
2615         }
2616     }
2617
2618     return hasVsite;
2619 }
2620
2621 //! Compute charge group centers of mass for molecule \p molt
2622 static void get_cgcm_mol(const gmx_moltype_t *molt,
2623                          const gmx_ffparams_t *ffparams,
2624                          int ePBC, t_graph *graph, const matrix box,
2625                          const rvec *x, rvec *xs, rvec *cg_cm)
2626 {
2627     int n, i;
2628
2629     if (ePBC != epbcNONE)
2630     {
2631         mk_mshift(nullptr, graph, ePBC, box, x);
2632
2633         shift_x(graph, box, x, xs);
2634         /* By doing an extra mk_mshift the molecules that are broken
2635          * because they were e.g. imported from another software
2636          * will be made whole again. Such are the healing powers
2637          * of GROMACS.
2638          */
2639         mk_mshift(nullptr, graph, ePBC, box, xs);
2640     }
2641     else
2642     {
2643         /* We copy the coordinates so the original coordinates remain
2644          * unchanged, just to be 100% sure that we do not affect
2645          * binary reproducibility of simulations.
2646          */
2647         n = molt->cgs.index[molt->cgs.nr];
2648         for (i = 0; i < n; i++)
2649         {
2650             copy_rvec(x[i], xs[i]);
2651         }
2652     }
2653
2654     if (moltypeHasVsite(*molt))
2655     {
2656         construct_vsites(nullptr, xs, 0.0, nullptr,
2657                          ffparams->iparams, molt->ilist,
2658                          epbcNONE, TRUE, nullptr, nullptr);
2659     }
2660
2661     calc_cgcm(nullptr, 0, molt->cgs.nr, &molt->cgs, xs, cg_cm);
2662 }
2663
2664 void dd_bonded_cg_distance(FILE *fplog,
2665                            const gmx_mtop_t *mtop,
2666                            const t_inputrec *ir,
2667                            const rvec *x, const matrix box,
2668                            gmx_bool bBCheck,
2669                            real *r_2b, real *r_mb)
2670 {
2671     gmx_bool           bExclRequired;
2672     int                at_offset;
2673     t_graph            graph;
2674     rvec              *xs, *cg_cm;
2675     bonded_distance_t  bd_2b = { 0, -1, -1, -1 };
2676     bonded_distance_t  bd_mb = { 0, -1, -1, -1 };
2677
2678     bExclRequired = inputrecExclForces(ir);
2679
2680     *r_2b     = 0;
2681     *r_mb     = 0;
2682     at_offset = 0;
2683     for (const gmx_molblock_t &molb : mtop->molblock)
2684     {
2685         const gmx_moltype_t &molt = mtop->moltype[molb.type];
2686         if (molt.cgs.nr == 1 || molb.nmol == 0)
2687         {
2688             at_offset += molb.nmol*molt.atoms.nr;
2689         }
2690         else
2691         {
2692             if (ir->ePBC != epbcNONE)
2693             {
2694                 mk_graph_ilist(nullptr, molt.ilist, 0, molt.atoms.nr, FALSE, FALSE,
2695                                &graph);
2696             }
2697
2698             std::vector<int> at2cg = make_at2cg(molt.cgs);
2699             snew(xs, molt.atoms.nr);
2700             snew(cg_cm, molt.cgs.nr);
2701             for (int mol = 0; mol < molb.nmol; mol++)
2702             {
2703                 get_cgcm_mol(&molt, &mtop->ffparams, ir->ePBC, &graph, box,
2704                              x+at_offset, xs, cg_cm);
2705
2706                 bonded_distance_t bd_mol_2b = { 0, -1, -1, -1 };
2707                 bonded_distance_t bd_mol_mb = { 0, -1, -1, -1 };
2708
2709                 bonded_cg_distance_mol(&molt, at2cg, bBCheck, bExclRequired, cg_cm,
2710                                        &bd_mol_2b, &bd_mol_mb);
2711
2712                 /* Process the mol data adding the atom index offset */
2713                 update_max_bonded_distance(bd_mol_2b.r2, bd_mol_2b.ftype,
2714                                            at_offset + bd_mol_2b.a1,
2715                                            at_offset + bd_mol_2b.a2,
2716                                            &bd_2b);
2717                 update_max_bonded_distance(bd_mol_mb.r2, bd_mol_mb.ftype,
2718                                            at_offset + bd_mol_mb.a1,
2719                                            at_offset + bd_mol_mb.a2,
2720                                            &bd_mb);
2721
2722                 at_offset += molt.atoms.nr;
2723             }
2724             sfree(cg_cm);
2725             sfree(xs);
2726             if (ir->ePBC != epbcNONE)
2727             {
2728                 done_graph(&graph);
2729             }
2730         }
2731     }
2732
2733     if (mtop->bIntermolecularInteractions)
2734     {
2735         if (ncg_mtop(mtop) < mtop->natoms)
2736         {
2737             gmx_fatal(FARGS, "The combination of intermolecular interactions, charge groups and domain decomposition is not supported. Use cutoff-scheme=Verlet (which removes the charge groups) or run without domain decomposition.");
2738         }
2739
2740         bonded_distance_intermol(mtop->intermolecular_ilist,
2741                                  bBCheck,
2742                                  x, ir->ePBC, box,
2743                                  &bd_2b, &bd_mb);
2744     }
2745
2746     *r_2b = sqrt(bd_2b.r2);
2747     *r_mb = sqrt(bd_mb.r2);
2748
2749     if (fplog && (*r_2b > 0 || *r_mb > 0))
2750     {
2751         fprintf(fplog,
2752                 "Initial maximum inter charge-group distances:\n");
2753         if (*r_2b > 0)
2754         {
2755             fprintf(fplog,
2756                     "    two-body bonded interactions: %5.3f nm, %s, atoms %d %d\n",
2757                     *r_2b, (bd_2b.ftype >= 0) ? interaction_function[bd_2b.ftype].longname : "Exclusion",
2758                     bd_2b.a1 + 1, bd_2b.a2 + 1);
2759         }
2760         if (*r_mb > 0)
2761         {
2762             fprintf(fplog,
2763                     "  multi-body bonded interactions: %5.3f nm, %s, atoms %d %d\n",
2764                     *r_mb, interaction_function[bd_mb.ftype].longname,
2765                     bd_mb.a1 + 1, bd_mb.a2 + 1);
2766         }
2767     }
2768 }