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