src/gromacs/domdec/domdec.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2005,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 #include "gmxpre.h"
  37
  38 #include "domdec.h"
  39
  40 #include "config.h"
  41
  42 #include <assert.h>
  43 #include <limits.h>
  44 #include <stdio.h>
  45 #include <stdlib.h>
  46 #include <string.h>
  47
  48 #include <cmath>
  49
  50 #include <algorithm>
  51
  52 #include "gromacs/compat/make_unique.h"
  53 #include "gromacs/domdec/collect.h"
  54 #include "gromacs/domdec/dlbtiming.h"
  55 #include "gromacs/domdec/domdec_network.h"
  56 #include "gromacs/domdec/ga2la.h"
  57 #include "gromacs/domdec/localatomsetmanager.h"
  58 #include "gromacs/ewald/pme.h"
  59 #include "gromacs/fileio/gmxfio.h"
  60 #include "gromacs/fileio/pdbio.h"
  61 #include "gromacs/gmxlib/chargegroup.h"
  62 #include "gromacs/gmxlib/network.h"
  63 #include "gromacs/gmxlib/nrnb.h"
  64 #include "gromacs/gpu_utils/gpu_utils.h"
  65 #include "gromacs/hardware/hw_info.h"
  66 #include "gromacs/imd/imd.h"
  67 #include "gromacs/listed-forces/manage-threading.h"
  68 #include "gromacs/math/functions.h"
  69 #include "gromacs/math/vec.h"
  70 #include "gromacs/math/vectypes.h"
  71 #include "gromacs/mdlib/constr.h"
  72 #include "gromacs/mdlib/constraintrange.h"
  73 #include "gromacs/mdlib/forcerec.h"
  74 #include "gromacs/mdlib/gmx_omp_nthreads.h"
  75 #include "gromacs/mdlib/lincs.h"
  76 #include "gromacs/mdlib/mdatoms.h"
  77 #include "gromacs/mdlib/mdrun.h"
  78 #include "gromacs/mdlib/mdsetup.h"
  79 #include "gromacs/mdlib/nb_verlet.h"
  80 #include "gromacs/mdlib/nbnxn_grid.h"
  81 #include "gromacs/mdlib/nsgrid.h"
  82 #include "gromacs/mdlib/vsite.h"
  83 #include "gromacs/mdtypes/commrec.h"
  84 #include "gromacs/mdtypes/df_history.h"
  85 #include "gromacs/mdtypes/forcerec.h"
  86 #include "gromacs/mdtypes/inputrec.h"
  87 #include "gromacs/mdtypes/md_enums.h"
  88 #include "gromacs/mdtypes/mdatom.h"
  89 #include "gromacs/mdtypes/nblist.h"
  90 #include "gromacs/mdtypes/state.h"
  91 #include "gromacs/pbcutil/ishift.h"
  92 #include "gromacs/pbcutil/pbc.h"
  93 #include "gromacs/pulling/pull.h"
  94 #include "gromacs/pulling/pull_rotation.h"
  95 #include "gromacs/timing/wallcycle.h"
  96 #include "gromacs/topology/block.h"
  97 #include "gromacs/topology/idef.h"
  98 #include "gromacs/topology/ifunc.h"
  99 #include "gromacs/topology/mtop_lookup.h"
 100 #include "gromacs/topology/mtop_util.h"
 101 #include "gromacs/topology/topology.h"
 102 #include "gromacs/utility/basedefinitions.h"
 103 #include "gromacs/utility/basenetwork.h"
 104 #include "gromacs/utility/cstringutil.h"
 105 #include "gromacs/utility/exceptions.h"
 106 #include "gromacs/utility/fatalerror.h"
 107 #include "gromacs/utility/gmxmpi.h"
 108 #include "gromacs/utility/qsort_threadsafe.h"
 109 #include "gromacs/utility/real.h"
 110 #include "gromacs/utility/smalloc.h"
 111 #include "gromacs/utility/strconvert.h"
 112 #include "gromacs/utility/stringutil.h"
 113
 114 #include "atomdistribution.h"
 115 #include "cellsizes.h"
 116 #include "distribute.h"
 117 #include "domdec_constraints.h"
 118 #include "domdec_internal.h"
 119 #include "domdec_vsite.h"
 120 #include "redistribute.h"
 121 #include "utility.h"
 122
 123 #define DD_NLOAD_MAX 9
 124
 125 static const char *edlbs_names[edlbsNR] = { "off", "auto", "locked", "on", "on" };
 126
 127 /* The size per charge group of the cggl_flag buffer in gmx_domdec_comm_t */
 128 #define DD_CGIBS 2
 129
 130 /* The flags for the cggl_flag buffer in gmx_domdec_comm_t */
 131 #define DD_FLAG_NRCG  65535
 132 #define DD_FLAG_FW(d) (1<<(16+(d)*2))
 133 #define DD_FLAG_BW(d) (1<<(16+(d)*2+1))
 134
 135 /* The DD zone order */
 136 static const ivec dd_zo[DD_MAXZONE] =
 137 {{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}, {0, 1, 1}, {0, 0, 1}, {1, 0, 1}, {1, 1, 1}};
 138
 139 /* The non-bonded zone-pair setup for domain decomposition
 140  * The first number is the i-zone, the second number the first j-zone seen by
 141  * this i-zone, the third number the last+1 j-zone seen by this i-zone.
 142  * As is, this is for 3D decomposition, where there are 4 i-zones.
 143  * With 2D decomposition use only the first 2 i-zones and a last+1 j-zone of 4.
 144  * With 1D decomposition use only the first i-zone and a last+1 j-zone of 2.
 145  */
 146 static const int
 147     ddNonbondedZonePairRanges[DD_MAXIZONE][3] = {{0, 0, 8},
 148                                                  {1, 3, 6},
 149                                                  {2, 5, 6},
 150                                                  {3, 5, 7}};
 151
 152 /* Turn on DLB when the load imbalance causes this amount of total loss.
 153  * There is a bit of overhead with DLB and it's difficult to achieve
 154  * a load imbalance of less than 2% with DLB.
 155  */
 156 #define DD_PERF_LOSS_DLB_ON  0.02
 157
 158 /* Warn about imbalance due to PP or PP/PME load imbalance at this loss */
 159 #define DD_PERF_LOSS_WARN    0.05
 160
 161
 162 /* We check if to turn on DLB at the first and every 100 DD partitionings.
 163  * With large imbalance DLB will turn on at the first step, so we can
 164  * make the interval so large that the MPI overhead of the check is negligible.
 165  */
 166 static const int c_checkTurnDlbOnInterval  = 100;
 167 /* We need to check if DLB results in worse performance and then turn it off.
 168  * We check this more often then for turning DLB on, because the DLB can scale
 169  * the domains very rapidly, so if unlucky the load imbalance can go up quickly
 170  * and furthermore, we are already synchronizing often with DLB, so
 171  * the overhead of the MPI Bcast is not that high.
 172  */
 173 static const int c_checkTurnDlbOffInterval =  20;
 174
 175 /* Forward declaration */
 176 static void dd_dlb_set_should_check_whether_to_turn_dlb_on(gmx_domdec_t *dd, gmx_bool bValue);
 177
 178
 179 /*
 180    #define dd_index(n,i) ((((i)[ZZ]*(n)[YY] + (i)[YY])*(n)[XX]) + (i)[XX])
 181
 182    static void index2xyz(ivec nc,int ind,ivec xyz)
 183    {
 184    xyz[XX] = ind % nc[XX];
 185    xyz[YY] = (ind / nc[XX]) % nc[YY];
 186    xyz[ZZ] = ind / (nc[YY]*nc[XX]);
 187    }
 188  */
 189
 190 static void ddindex2xyz(const ivec nc, int ind, ivec xyz)
 191 {
 192     xyz[XX] = ind / (nc[YY]*nc[ZZ]);
 193     xyz[YY] = (ind / nc[ZZ]) % nc[YY];
 194     xyz[ZZ] = ind % nc[ZZ];
 195 }
 196
 197 static int ddcoord2ddnodeid(gmx_domdec_t *dd, ivec c)
 198 {
 199     int ddindex;
 200     int ddnodeid = -1;
 201
 202     ddindex = dd_index(dd->nc, c);
 203     if (dd->comm->bCartesianPP_PME)
 204     {
 205         ddnodeid = dd->comm->ddindex2ddnodeid[ddindex];
 206     }
 207     else if (dd->comm->bCartesianPP)
 208     {
 209 #if GMX_MPI
 210         MPI_Cart_rank(dd->mpi_comm_all, c, &ddnodeid);
 211 #endif
 212     }
 213     else
 214     {
 215         ddnodeid = ddindex;
 216     }
 217
 218     return ddnodeid;
 219 }
 220
 221 static gmx_bool dynamic_dd_box(const gmx_ddbox_t *ddbox, const t_inputrec *ir)
 222 {
 223     return (ddbox->nboundeddim < DIM || inputrecDynamicBox(ir));
 224 }
 225
 226 int ddglatnr(const gmx_domdec_t *dd, int i)
 227 {
 228     int atnr;
 229
 230     if (dd == nullptr)
 231     {
 232         atnr = i + 1;
 233     }
 234     else
 235     {
 236         if (i >= dd->comm->atomRanges.numAtomsTotal())
 237         {
 238             gmx_fatal(FARGS, "glatnr called with %d, which is larger than the local number of atoms (%d)", i, dd->comm->atomRanges.numAtomsTotal());
 239         }
 240         atnr = dd->globalAtomIndices[i] + 1;
 241     }
 242
 243     return atnr;
 244 }
 245
 246 t_block *dd_charge_groups_global(gmx_domdec_t *dd)
 247 {
 248     return &dd->comm->cgs_gl;
 249 }
 250
 251 void dd_store_state(gmx_domdec_t *dd, t_state *state)
 252 {
 253     int i;
 254
 255     if (state->ddp_count != dd->ddp_count)
 256     {
 257         gmx_incons("The MD state does not match the domain decomposition state");
 258     }
 259
 260     state->cg_gl.resize(dd->ncg_home);
 261     for (i = 0; i < dd->ncg_home; i++)
 262     {
 263         state->cg_gl[i] = dd->globalAtomGroupIndices[i];
 264     }
 265
 266     state->ddp_count_cg_gl = dd->ddp_count;
 267 }
 268
 269 gmx_domdec_zones_t *domdec_zones(gmx_domdec_t *dd)
 270 {
 271     return &dd->comm->zones;
 272 }
 273
 274 void dd_get_ns_ranges(const gmx_domdec_t *dd, int icg,
 275                       int *jcg0, int *jcg1, ivec shift0, ivec shift1)
 276 {
 277     gmx_domdec_zones_t *zones;
 278     int                 izone, d, dim;
 279
 280     zones = &dd->comm->zones;
 281
 282     izone = 0;
 283     while (icg >= zones->izone[izone].cg1)
 284     {
 285         izone++;
 286     }
 287
 288     if (izone == 0)
 289     {
 290         *jcg0 = icg;
 291     }
 292     else if (izone < zones->nizone)
 293     {
 294         *jcg0 = zones->izone[izone].jcg0;
 295     }
 296     else
 297     {
 298         gmx_fatal(FARGS, "DD icg %d out of range: izone (%d) >= nizone (%d)",
 299                   icg, izone, zones->nizone);
 300     }
 301
 302     *jcg1 = zones->izone[izone].jcg1;
 303
 304     for (d = 0; d < dd->ndim; d++)
 305     {
 306         dim         = dd->dim[d];
 307         shift0[dim] = zones->izone[izone].shift0[dim];
 308         shift1[dim] = zones->izone[izone].shift1[dim];
 309         if (dd->comm->tric_dir[dim] || (isDlbOn(dd->comm) && d > 0))
 310         {
 311             /* A conservative approach, this can be optimized */
 312             shift0[dim] -= 1;
 313             shift1[dim] += 1;
 314         }
 315     }
 316 }
 317
 318 int dd_numHomeAtoms(const gmx_domdec_t &dd)
 319 {
 320     return dd.comm->atomRanges.numHomeAtoms();
 321 }
 322
 323 int dd_natoms_mdatoms(const gmx_domdec_t *dd)
 324 {
 325     /* We currently set mdatoms entries for all atoms:
 326      * local + non-local + communicated for vsite + constraints
 327      */
 328
 329     return dd->comm->atomRanges.numAtomsTotal();
 330 }
 331
 332 int dd_natoms_vsite(const gmx_domdec_t *dd)
 333 {
 334     return dd->comm->atomRanges.end(DDAtomRanges::Type::Vsites);
 335 }
 336
 337 void dd_get_constraint_range(const gmx_domdec_t *dd, int *at_start, int *at_end)
 338 {
 339     *at_start = dd->comm->atomRanges.start(DDAtomRanges::Type::Constraints);
 340     *at_end   = dd->comm->atomRanges.end(DDAtomRanges::Type::Constraints);
 341 }
 342
 343 void dd_move_x(gmx_domdec_t             *dd,
 344                matrix                    box,
 345                gmx::ArrayRef<gmx::RVec>  x,
 346                gmx_wallcycle            *wcycle)
 347 {
 348     wallcycle_start(wcycle, ewcMOVEX);
 349
 350     int                    nzone, nat_tot;
 351     gmx_domdec_comm_t     *comm;
 352     gmx_domdec_comm_dim_t *cd;
 353     rvec                   shift = {0, 0, 0};
 354     gmx_bool               bPBC, bScrew;
 355
 356     comm = dd->comm;
 357
 358     const gmx::RangePartitioning &atomGrouping = dd->atomGrouping();
 359
 360     nzone   = 1;
 361     nat_tot = comm->atomRanges.numHomeAtoms();
 362     for (int d = 0; d < dd->ndim; d++)
 363     {
 364         bPBC   = (dd->ci[dd->dim[d]] == 0);
 365         bScrew = (bPBC && dd->bScrewPBC && dd->dim[d] == XX);
 366         if (bPBC)
 367         {
 368             copy_rvec(box[dd->dim[d]], shift);
 369         }
 370         cd = &comm->cd[d];
 371         for (const gmx_domdec_ind_t &ind : cd->ind)
 372         {
 373             DDBufferAccess<gmx::RVec>  sendBufferAccess(comm->rvecBuffer, ind.nsend[nzone + 1]);
 374             gmx::ArrayRef<gmx::RVec>  &sendBuffer = sendBufferAccess.buffer;
 375             int                        n          = 0;
 376             if (!bPBC)
 377             {
 378                 for (int g : ind.index)
 379                 {
 380                     for (int j : atomGrouping.block(g))
 381                     {
 382                         sendBuffer[n] = x[j];
 383                         n++;
 384                     }
 385                 }
 386             }
 387             else if (!bScrew)
 388             {
 389                 for (int g : ind.index)
 390                 {
 391                     for (int j : atomGrouping.block(g))
 392                     {
 393                         /* We need to shift the coordinates */
 394                         for (int d = 0; d < DIM; d++)
 395                         {
 396                             sendBuffer[n][d] = x[j][d] + shift[d];
 397                         }
 398                         n++;
 399                     }
 400                 }
 401             }
 402             else
 403             {
 404                 for (int g : ind.index)
 405                 {
 406                     for (int j : atomGrouping.block(g))
 407                     {
 408                         /* Shift x */
 409                         sendBuffer[n][XX] = x[j][XX] + shift[XX];
 410                         /* Rotate y and z.
 411                          * This operation requires a special shift force
 412                          * treatment, which is performed in calc_vir.
 413                          */
 414                         sendBuffer[n][YY] = box[YY][YY] - x[j][YY];
 415                         sendBuffer[n][ZZ] = box[ZZ][ZZ] - x[j][ZZ];
 416                         n++;
 417                     }
 418                 }
 419             }
 420
 421             DDBufferAccess<gmx::RVec>  receiveBufferAccess(comm->rvecBuffer2, cd->receiveInPlace ? 0 : ind.nrecv[nzone + 1]);
 422
 423             gmx::ArrayRef<gmx::RVec>   receiveBuffer;
 424             if (cd->receiveInPlace)
 425             {
 426                 receiveBuffer = gmx::arrayRefFromArray(x.data() + nat_tot, ind.nrecv[nzone + 1]);
 427             }
 428             else
 429             {
 430                 receiveBuffer = receiveBufferAccess.buffer;
 431             }
 432             /* Send and receive the coordinates */
 433             ddSendrecv(dd, d, dddirBackward,
 434                        sendBuffer, receiveBuffer);
 435
 436             if (!cd->receiveInPlace)
 437             {
 438                 int j = 0;
 439                 for (int zone = 0; zone < nzone; zone++)
 440                 {
 441                     for (int i = ind.cell2at0[zone]; i < ind.cell2at1[zone]; i++)
 442                     {
 443                         x[i] = receiveBuffer[j++];
 444                     }
 445                 }
 446             }
 447             nat_tot += ind.nrecv[nzone+1];
 448         }
 449         nzone += nzone;
 450     }
 451
 452     wallcycle_stop(wcycle, ewcMOVEX);
 453 }
 454
 455 void dd_move_f(gmx_domdec_t             *dd,
 456                gmx::ArrayRef<gmx::RVec>  f,
 457                rvec                     *fshift,
 458                gmx_wallcycle            *wcycle)
 459 {
 460     wallcycle_start(wcycle, ewcMOVEF);
 461
 462     int                    nzone, nat_tot;
 463     gmx_domdec_comm_t     *comm;
 464     gmx_domdec_comm_dim_t *cd;
 465     ivec                   vis;
 466     int                    is;
 467     gmx_bool               bShiftForcesNeedPbc, bScrew;
 468
 469     comm = dd->comm;
 470
 471     const gmx::RangePartitioning &atomGrouping = dd->atomGrouping();
 472
 473     nzone   = comm->zones.n/2;
 474     nat_tot = comm->atomRanges.end(DDAtomRanges::Type::Zones);
 475     for (int d = dd->ndim-1; d >= 0; d--)
 476     {
 477         /* Only forces in domains near the PBC boundaries need to
 478            consider PBC in the treatment of fshift */
 479         bShiftForcesNeedPbc   = (dd->ci[dd->dim[d]] == 0);
 480         bScrew                = (bShiftForcesNeedPbc && dd->bScrewPBC && dd->dim[d] == XX);
 481         if (fshift == nullptr && !bScrew)
 482         {
 483             bShiftForcesNeedPbc = FALSE;
 484         }
 485         /* Determine which shift vector we need */
 486         clear_ivec(vis);
 487         vis[dd->dim[d]] = 1;
 488         is              = IVEC2IS(vis);
 489
 490         cd = &comm->cd[d];
 491         for (int p = cd->numPulses() - 1; p >= 0; p--)
 492         {
 493             const gmx_domdec_ind_t    &ind  = cd->ind[p];
 494             DDBufferAccess<gmx::RVec>  receiveBufferAccess(comm->rvecBuffer, ind.nsend[nzone + 1]);
 495             gmx::ArrayRef<gmx::RVec>  &receiveBuffer = receiveBufferAccess.buffer;
 496
 497             nat_tot                        -= ind.nrecv[nzone+1];
 498
 499             DDBufferAccess<gmx::RVec>  sendBufferAccess(comm->rvecBuffer2, cd->receiveInPlace ? 0 : ind.nrecv[nzone + 1]);
 500
 501             gmx::ArrayRef<gmx::RVec>   sendBuffer;
 502             if (cd->receiveInPlace)
 503             {
 504                 sendBuffer = gmx::arrayRefFromArray(f.data() + nat_tot, ind.nrecv[nzone + 1]);
 505             }
 506             else
 507             {
 508                 sendBuffer = sendBufferAccess.buffer;
 509                 int j = 0;
 510                 for (int zone = 0; zone < nzone; zone++)
 511                 {
 512                     for (int i = ind.cell2at0[zone]; i < ind.cell2at1[zone]; i++)
 513                     {
 514                         sendBuffer[j++] = f[i];
 515                     }
 516                 }
 517             }
 518             /* Communicate the forces */
 519             ddSendrecv(dd, d, dddirForward,
 520                        sendBuffer, receiveBuffer);
 521             /* Add the received forces */
 522             int n = 0;
 523             if (!bShiftForcesNeedPbc)
 524             {
 525                 for (int g : ind.index)
 526                 {
 527                     for (int j : atomGrouping.block(g))
 528                     {
 529                         for (int d = 0; d < DIM; d++)
 530                         {
 531                             f[j][d] += receiveBuffer[n][d];
 532                         }
 533                         n++;
 534                     }
 535                 }
 536             }
 537             else if (!bScrew)
 538             {
 539                 /* fshift should always be defined if this function is
 540                  * called when bShiftForcesNeedPbc is true */
 541                 assert(nullptr != fshift);
 542                 for (int g : ind.index)
 543                 {
 544                     for (int j : atomGrouping.block(g))
 545                     {
 546                         for (int d = 0; d < DIM; d++)
 547                         {
 548                             f[j][d] += receiveBuffer[n][d];
 549                         }
 550                         /* Add this force to the shift force */
 551                         for (int d = 0; d < DIM; d++)
 552                         {
 553                             fshift[is][d] += receiveBuffer[n][d];
 554                         }
 555                         n++;
 556                     }
 557                 }
 558             }
 559             else
 560             {
 561                 for (int g : ind.index)
 562                 {
 563                     for (int j : atomGrouping.block(g))
 564                     {
 565                         /* Rotate the force */
 566                         f[j][XX] += receiveBuffer[n][XX];
 567                         f[j][YY] -= receiveBuffer[n][YY];
 568                         f[j][ZZ] -= receiveBuffer[n][ZZ];
 569                         if (fshift)
 570                         {
 571                             /* Add this force to the shift force */
 572                             for (int d = 0; d < DIM; d++)
 573                             {
 574                                 fshift[is][d] += receiveBuffer[n][d];
 575                             }
 576                         }
 577                         n++;
 578                     }
 579                 }
 580             }
 581         }
 582         nzone /= 2;
 583     }
 584     wallcycle_stop(wcycle, ewcMOVEF);
 585 }
 586
 587 /* Convenience function for extracting a real buffer from an rvec buffer
 588  *
 589  * To reduce the number of temporary communication buffers and avoid
 590  * cache polution, we reuse gmx::RVec buffers for storing reals.
 591  * This functions return a real buffer reference with the same number
 592  * of elements as the gmx::RVec buffer (so 1/3 of the size in bytes).
 593  */
 594 static gmx::ArrayRef<real>
 595 realArrayRefFromRvecArrayRef(gmx::ArrayRef<gmx::RVec> arrayRef)
 596 {
 597     return gmx::arrayRefFromArray(as_rvec_array(arrayRef.data())[0],
 598                                   arrayRef.size());
 599 }
 600
 601 void dd_atom_spread_real(gmx_domdec_t *dd, real v[])
 602 {
 603     int                    nzone, nat_tot;
 604     gmx_domdec_comm_t     *comm;
 605     gmx_domdec_comm_dim_t *cd;
 606
 607     comm = dd->comm;
 608
 609     const gmx::RangePartitioning &atomGrouping = dd->atomGrouping();
 610
 611     nzone   = 1;
 612     nat_tot = comm->atomRanges.numHomeAtoms();
 613     for (int d = 0; d < dd->ndim; d++)
 614     {
 615         cd = &comm->cd[d];
 616         for (const gmx_domdec_ind_t &ind : cd->ind)
 617         {
 618             /* Note: We provision for RVec instead of real, so a factor of 3
 619              * more than needed. The buffer actually already has this size
 620              * and we pass a plain pointer below, so this does not matter.
 621              */
 622             DDBufferAccess<gmx::RVec> sendBufferAccess(comm->rvecBuffer, ind.nsend[nzone + 1]);
 623             gmx::ArrayRef<real>       sendBuffer = realArrayRefFromRvecArrayRef(sendBufferAccess.buffer);
 624             int                       n          = 0;
 625             for (int g : ind.index)
 626             {
 627                 for (int j : atomGrouping.block(g))
 628                 {
 629                     sendBuffer[n++] = v[j];
 630                 }
 631             }
 632
 633             DDBufferAccess<gmx::RVec> receiveBufferAccess(comm->rvecBuffer2, cd->receiveInPlace ? 0 : ind.nrecv[nzone + 1]);
 634
 635             gmx::ArrayRef<real>       receiveBuffer;
 636             if (cd->receiveInPlace)
 637             {
 638                 receiveBuffer = gmx::arrayRefFromArray(v + nat_tot, ind.nrecv[nzone + 1]);
 639             }
 640             else
 641             {
 642                 receiveBuffer = realArrayRefFromRvecArrayRef(receiveBufferAccess.buffer);
 643             }
 644             /* Send and receive the data */
 645             ddSendrecv(dd, d, dddirBackward,
 646                        sendBuffer, receiveBuffer);
 647             if (!cd->receiveInPlace)
 648             {
 649                 int j = 0;
 650                 for (int zone = 0; zone < nzone; zone++)
 651                 {
 652                     for (int i = ind.cell2at0[zone]; i < ind.cell2at1[zone]; i++)
 653                     {
 654                         v[i] = receiveBuffer[j++];
 655                     }
 656                 }
 657             }
 658             nat_tot += ind.nrecv[nzone+1];
 659         }
 660         nzone += nzone;
 661     }
 662 }
 663
 664 void dd_atom_sum_real(gmx_domdec_t *dd, real v[])
 665 {
 666     int                    nzone, nat_tot;
 667     gmx_domdec_comm_t     *comm;
 668     gmx_domdec_comm_dim_t *cd;
 669
 670     comm = dd->comm;
 671
 672     const gmx::RangePartitioning &atomGrouping = dd->atomGrouping();
 673
 674     nzone   = comm->zones.n/2;
 675     nat_tot = comm->atomRanges.end(DDAtomRanges::Type::Zones);
 676     for (int d = dd->ndim-1; d >= 0; d--)
 677     {
 678         cd = &comm->cd[d];
 679         for (int p = cd->numPulses() - 1; p >= 0; p--)
 680         {
 681             const gmx_domdec_ind_t &ind = cd->ind[p];
 682
 683             /* Note: We provision for RVec instead of real, so a factor of 3
 684              * more than needed. The buffer actually already has this size
 685              * and we typecast, so this works as intended.
 686              */
 687             DDBufferAccess<gmx::RVec> receiveBufferAccess(comm->rvecBuffer, ind.nsend[nzone + 1]);
 688             gmx::ArrayRef<real>       receiveBuffer = realArrayRefFromRvecArrayRef(receiveBufferAccess.buffer);
 689             nat_tot -= ind.nrecv[nzone + 1];
 690
 691             DDBufferAccess<gmx::RVec> sendBufferAccess(comm->rvecBuffer2, cd->receiveInPlace ? 0 : ind.nrecv[nzone + 1]);
 692
 693             gmx::ArrayRef<real>       sendBuffer;
 694             if (cd->receiveInPlace)
 695             {
 696                 sendBuffer = gmx::arrayRefFromArray(v + nat_tot, ind.nrecv[nzone + 1]);
 697             }
 698             else
 699             {
 700                 sendBuffer = realArrayRefFromRvecArrayRef(sendBufferAccess.buffer);
 701                 int j = 0;
 702                 for (int zone = 0; zone < nzone; zone++)
 703                 {
 704                     for (int i = ind.cell2at0[zone]; i < ind.cell2at1[zone]; i++)
 705                     {
 706                         sendBuffer[j++] = v[i];
 707                     }
 708                 }
 709             }
 710             /* Communicate the forces */
 711             ddSendrecv(dd, d, dddirForward,
 712                        sendBuffer, receiveBuffer);
 713             /* Add the received forces */
 714             int n = 0;
 715             for (int g : ind.index)
 716             {
 717                 for (int j : atomGrouping.block(g))
 718                 {
 719                     v[j] += receiveBuffer[n];
 720                     n++;
 721                 }
 722             }
 723         }
 724         nzone /= 2;
 725     }
 726 }
 727
 728 static void print_ddzone(FILE *fp, int d, int i, int j, gmx_ddzone_t *zone)
 729 {
 730     fprintf(fp, "zone d0 %d d1 %d d2 %d  min0 %6.3f max1 %6.3f mch0 %6.3f mch1 %6.3f p1_0 %6.3f p1_1 %6.3f\n",
 731             d, i, j,
 732             zone->min0, zone->max1,
 733             zone->mch0, zone->mch0,
 734             zone->p1_0, zone->p1_1);
 735 }
 736
 737 /* Using the home grid size as input in cell_ns_x0 and cell_ns_x1
 738  * takes the extremes over all home and remote zones in the halo
 739  * and returns the results in cell_ns_x0 and cell_ns_x1.
 740  * Note: only used with the group cut-off scheme.
 741  */
 742 static void dd_move_cellx(gmx_domdec_t      *dd,
 743                           const gmx_ddbox_t *ddbox,
 744                           rvec               cell_ns_x0,
 745                           rvec               cell_ns_x1)
 746 {
 747     constexpr int      c_ddZoneCommMaxNumZones = 5;
 748     gmx_ddzone_t       buf_s[c_ddZoneCommMaxNumZones];
 749     gmx_ddzone_t       buf_r[c_ddZoneCommMaxNumZones];
 750     gmx_ddzone_t       buf_e[c_ddZoneCommMaxNumZones];
 751     gmx_domdec_comm_t *comm = dd->comm;
 752
 753     rvec               extr_s[2];
 754     rvec               extr_r[2];
 755     for (int d = 1; d < dd->ndim; d++)
 756     {
 757         int           dim = dd->dim[d];
 758         gmx_ddzone_t &zp  = (d == 1) ? comm->zone_d1[0] : comm->zone_d2[0][0];
 759
 760         /* Copy the base sizes of the home zone */
 761         zp.min0 = cell_ns_x0[dim];
 762         zp.max1 = cell_ns_x1[dim];
 763         zp.min1 = cell_ns_x1[dim];
 764         zp.mch0 = cell_ns_x0[dim];
 765         zp.mch1 = cell_ns_x1[dim];
 766         zp.p1_0 = cell_ns_x0[dim];
 767         zp.p1_1 = cell_ns_x1[dim];
 768     }
 769
 770     gmx::ArrayRef<DDCellsizesWithDlb> cellsizes = comm->cellsizesWithDlb;
 771
 772     /* Loop backward over the dimensions and aggregate the extremes
 773      * of the cell sizes.
 774      */
 775     for (int d = dd->ndim - 2; d >= 0; d--)
 776     {
 777         const int  dim      = dd->dim[d];
 778         const bool applyPbc = (dim < ddbox->npbcdim);
 779
 780         /* Use an rvec to store two reals */
 781         extr_s[d][0] = cellsizes[d + 1].fracLower;
 782         extr_s[d][1] = cellsizes[d + 1].fracUpper;
 783         extr_s[d][2] = cellsizes[d + 1].fracUpper;
 784
 785         int pos = 0;
 786         GMX_ASSERT(pos < c_ddZoneCommMaxNumZones, "The buffers should be sufficiently large");
 787         /* Store the extremes in the backward sending buffer,
 788          * so they get updated separately from the forward communication.
 789          */
 790         for (int d1 = d; d1 < dd->ndim-1; d1++)
 791         {
 792             /* We invert the order to be able to use the same loop for buf_e */
 793             buf_s[pos].min0 = extr_s[d1][1];
 794             buf_s[pos].max1 = extr_s[d1][0];
 795             buf_s[pos].min1 = extr_s[d1][2];
 796             buf_s[pos].mch0 = 0;
 797             buf_s[pos].mch1 = 0;
 798             /* Store the cell corner of the dimension we communicate along */
 799             buf_s[pos].p1_0 = comm->cell_x0[dim];
 800             buf_s[pos].p1_1 = 0;
 801             pos++;
 802         }
 803
 804         buf_s[pos] = (dd->ndim == 2) ? comm->zone_d1[0] : comm->zone_d2[0][0];
 805         pos++;
 806
 807         if (dd->ndim == 3 && d == 0)
 808         {
 809             buf_s[pos] = comm->zone_d2[0][1];
 810             pos++;
 811             buf_s[pos] = comm->zone_d1[0];
 812             pos++;
 813         }
 814
 815         /* We only need to communicate the extremes
 816          * in the forward direction
 817          */
 818         int numPulses = comm->cd[d].numPulses();
 819         int numPulsesMin;
 820         if (applyPbc)
 821         {
 822             /* Take the minimum to avoid double communication */
 823             numPulsesMin = std::min(numPulses, dd->nc[dim] - 1 - numPulses);
 824         }
 825         else
 826         {
 827             /* Without PBC we should really not communicate over
 828              * the boundaries, but implementing that complicates
 829              * the communication setup and therefore we simply
 830              * do all communication, but ignore some data.
 831              */
 832             numPulsesMin = numPulses;
 833         }
 834         for (int pulse = 0; pulse < numPulsesMin; pulse++)
 835         {
 836             /* Communicate the extremes forward */
 837             bool receiveValidData = (applyPbc || dd->ci[dim] > 0);
 838
 839             int  numElements      = dd->ndim - d - 1;
 840             ddSendrecv(dd, d, dddirForward,
 841                        extr_s + d, numElements,
 842                        extr_r + d, numElements);
 843
 844             if (receiveValidData)
 845             {
 846                 for (int d1 = d; d1 < dd->ndim - 1; d1++)
 847                 {
 848                     extr_s[d1][0] = std::max(extr_s[d1][0], extr_r[d1][0]);
 849                     extr_s[d1][1] = std::min(extr_s[d1][1], extr_r[d1][1]);
 850                     extr_s[d1][2] = std::min(extr_s[d1][2], extr_r[d1][2]);
 851                 }
 852             }
 853         }
 854
 855         const int numElementsInBuffer = pos;
 856         for (int pulse = 0; pulse < numPulses; pulse++)
 857         {
 858             /* Communicate all the zone information backward */
 859             bool receiveValidData = (applyPbc || dd->ci[dim] < dd->nc[dim] - 1);
 860
 861             static_assert(sizeof(gmx_ddzone_t) == c_ddzoneNumReals*sizeof(real), "Here we expect gmx_ddzone_t to consist of c_ddzoneNumReals reals (only)");
 862
 863             int numReals = numElementsInBuffer*c_ddzoneNumReals;
 864             ddSendrecv(dd, d, dddirBackward,
 865                        gmx::arrayRefFromArray(&buf_s[0].min0, numReals),
 866                        gmx::arrayRefFromArray(&buf_r[0].min0, numReals));
 867
 868             rvec dh = { 0 };
 869             if (pulse > 0)
 870             {
 871                 for (int d1 = d + 1; d1 < dd->ndim; d1++)
 872                 {
 873                     /* Determine the decrease of maximum required
 874                      * communication height along d1 due to the distance along d,
 875                      * this avoids a lot of useless atom communication.
 876                      */
 877                     real dist_d = comm->cell_x1[dim] - buf_r[0].p1_0;
 878
 879                     int  c;
 880                     if (ddbox->tric_dir[dim])
 881                     {
 882                         /* c is the off-diagonal coupling between the cell planes
 883                          * along directions d and d1.
 884                          */
 885                         c = ddbox->v[dim][dd->dim[d1]][dim];
 886                     }
 887                     else
 888                     {
 889                         c = 0;
 890                     }
 891                     real det = (1 + c*c)*comm->cutoff*comm->cutoff - dist_d*dist_d;
 892                     if (det > 0)
 893                     {
 894                         dh[d1] = comm->cutoff - (c*dist_d + std::sqrt(det))/(1 + c*c);
 895                     }
 896                     else
 897                     {
 898                         /* A negative value signals out of range */
 899                         dh[d1] = -1;
 900                     }
 901                 }
 902             }
 903
 904             /* Accumulate the extremes over all pulses */
 905             for (int i = 0; i < numElementsInBuffer; i++)
 906             {
 907                 if (pulse == 0)
 908                 {
 909                     buf_e[i] = buf_r[i];
 910                 }
 911                 else
 912                 {
 913                     if (receiveValidData)
 914                     {
 915                         buf_e[i].min0 = std::min(buf_e[i].min0, buf_r[i].min0);
 916                         buf_e[i].max1 = std::max(buf_e[i].max1, buf_r[i].max1);
 917                         buf_e[i].min1 = std::min(buf_e[i].min1, buf_r[i].min1);
 918                     }
 919
 920                     int d1;
 921                     if (dd->ndim == 3 && d == 0 && i == numElementsInBuffer - 1)
 922                     {
 923                         d1 = 1;
 924                     }
 925                     else
 926                     {
 927                         d1 = d + 1;
 928                     }
 929                     if (receiveValidData && dh[d1] >= 0)
 930                     {
 931                         buf_e[i].mch0 = std::max(buf_e[i].mch0, buf_r[i].mch0-dh[d1]);
 932                         buf_e[i].mch1 = std::max(buf_e[i].mch1, buf_r[i].mch1-dh[d1]);
 933                     }
 934                 }
 935                 /* Copy the received buffer to the send buffer,
 936                  * to pass the data through with the next pulse.
 937                  */
 938                 buf_s[i] = buf_r[i];
 939             }
 940             if (((applyPbc || dd->ci[dim] + numPulses < dd->nc[dim]) && pulse == numPulses - 1) ||
 941                 (!applyPbc && dd->ci[dim] + 1 + pulse == dd->nc[dim] - 1))
 942             {
 943                 /* Store the extremes */
 944                 int pos = 0;
 945
 946                 for (int d1 = d; d1 < dd->ndim-1; d1++)
 947                 {
 948                     extr_s[d1][1] = std::min(extr_s[d1][1], buf_e[pos].min0);
 949                     extr_s[d1][0] = std::max(extr_s[d1][0], buf_e[pos].max1);
 950                     extr_s[d1][2] = std::min(extr_s[d1][2], buf_e[pos].min1);
 951                     pos++;
 952                 }
 953
 954                 if (d == 1 || (d == 0 && dd->ndim == 3))
 955                 {
 956                     for (int i = d; i < 2; i++)
 957                     {
 958                         comm->zone_d2[1-d][i] = buf_e[pos];
 959                         pos++;
 960                     }
 961                 }
 962                 if (d == 0)
 963                 {
 964                     comm->zone_d1[1] = buf_e[pos];
 965                     pos++;
 966                 }
 967             }
 968         }
 969     }
 970
 971     if (dd->ndim >= 2)
 972     {
 973         int dim = dd->dim[1];
 974         for (int i = 0; i < 2; i++)
 975         {
 976             if (debug)
 977             {
 978                 print_ddzone(debug, 1, i, 0, &comm->zone_d1[i]);
 979             }
 980             cell_ns_x0[dim] = std::min(cell_ns_x0[dim], comm->zone_d1[i].min0);
 981             cell_ns_x1[dim] = std::max(cell_ns_x1[dim], comm->zone_d1[i].max1);
 982         }
 983     }
 984     if (dd->ndim >= 3)
 985     {
 986         int dim = dd->dim[2];
 987         for (int i = 0; i < 2; i++)
 988         {
 989             for (int j = 0; j < 2; j++)
 990             {
 991                 if (debug)
 992                 {
 993                     print_ddzone(debug, 2, i, j, &comm->zone_d2[i][j]);
 994                 }
 995                 cell_ns_x0[dim] = std::min(cell_ns_x0[dim], comm->zone_d2[i][j].min0);
 996                 cell_ns_x1[dim] = std::max(cell_ns_x1[dim], comm->zone_d2[i][j].max1);
 997             }
 998         }
 999     }
1000     for (int d = 1; d < dd->ndim; d++)
1001     {
1002         cellsizes[d].fracLowerMax = extr_s[d-1][0];
1003         cellsizes[d].fracUpperMin = extr_s[d-1][1];
1004         if (debug)
1005         {
1006             fprintf(debug, "Cell fraction d %d, max0 %f, min1 %f\n",
1007                     d, cellsizes[d].fracLowerMax, cellsizes[d].fracUpperMin);
1008         }
1009     }
1010 }
1011
1012 static void write_dd_grid_pdb(const char *fn, int64_t step,
1013                               gmx_domdec_t *dd, matrix box, gmx_ddbox_t *ddbox)
1014 {
1015     rvec   grid_s[2], *grid_r = nullptr, cx, r;
1016     char   fname[STRLEN], buf[22];
1017     FILE  *out;
1018     int    a, i, d, z, y, x;
1019     matrix tric;
1020     real   vol;
1021
1022     copy_rvec(dd->comm->cell_x0, grid_s[0]);
1023     copy_rvec(dd->comm->cell_x1, grid_s[1]);
1024
1025     if (DDMASTER(dd))
1026     {
1027         snew(grid_r, 2*dd->nnodes);
1028     }
1029
1030     dd_gather(dd, 2*sizeof(rvec), grid_s, DDMASTER(dd) ? grid_r : nullptr);
1031
1032     if (DDMASTER(dd))
1033     {
1034         for (d = 0; d < DIM; d++)
1035         {
1036             for (i = 0; i < DIM; i++)
1037             {
1038                 if (d == i)
1039                 {
1040                     tric[d][i] = 1;
1041                 }
1042                 else
1043                 {
1044                     if (d < ddbox->npbcdim && dd->nc[d] > 1)
1045                     {
1046                         tric[d][i] = box[i][d]/box[i][i];
1047                     }
1048                     else
1049                     {
1050                         tric[d][i] = 0;
1051                     }
1052                 }
1053             }
1054         }
1055         sprintf(fname, "%s_%s.pdb", fn, gmx_step_str(step, buf));
1056         out = gmx_fio_fopen(fname, "w");
1057         gmx_write_pdb_box(out, dd->bScrewPBC ? epbcSCREW : epbcXYZ, box);
1058         a = 1;
1059         for (i = 0; i < dd->nnodes; i++)
1060         {
1061             vol = dd->nnodes/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
1062             for (d = 0; d < DIM; d++)
1063             {
1064                 vol *= grid_r[i*2+1][d] - grid_r[i*2][d];
1065             }
1066             for (z = 0; z < 2; z++)
1067             {
1068                 for (y = 0; y < 2; y++)
1069                 {
1070                     for (x = 0; x < 2; x++)
1071                     {
1072                         cx[XX] = grid_r[i*2+x][XX];
1073                         cx[YY] = grid_r[i*2+y][YY];
1074                         cx[ZZ] = grid_r[i*2+z][ZZ];
1075                         mvmul(tric, cx, r);
1076                         gmx_fprintf_pdb_atomline(out, epdbATOM, a++, "CA", ' ', "GLY", ' ', i+1, ' ',
1077                                                  10*r[XX], 10*r[YY], 10*r[ZZ], 1.0, vol, "");
1078                     }
1079                 }
1080             }
1081             for (d = 0; d < DIM; d++)
1082             {
1083                 for (x = 0; x < 4; x++)
1084                 {
1085                     switch (d)
1086                     {
1087                         case 0: y = 1 + i*8 + 2*x; break;
1088                         case 1: y = 1 + i*8 + 2*x - (x % 2); break;
1089                         case 2: y = 1 + i*8 + x; break;
1090                     }
1091                     fprintf(out, "%6s%5d%5d\n", "CONECT", y, y+(1<<d));
1092                 }
1093             }
1094         }
1095         gmx_fio_fclose(out);
1096         sfree(grid_r);
1097     }
1098 }
1099
1100 void write_dd_pdb(const char *fn, int64_t step, const char *title,
1101                   const gmx_mtop_t *mtop, const t_commrec *cr,
1102                   int natoms, const rvec x[], const matrix box)
1103 {
1104     char          fname[STRLEN], buf[22];
1105     FILE         *out;
1106     int           resnr;
1107     const char   *atomname, *resname;
1108     gmx_domdec_t *dd;
1109
1110     dd = cr->dd;
1111     if (natoms == -1)
1112     {
1113         natoms = dd->comm->atomRanges.end(DDAtomRanges::Type::Vsites);
1114     }
1115
1116     sprintf(fname, "%s_%s_n%d.pdb", fn, gmx_step_str(step, buf), cr->sim_nodeid);
1117
1118     out = gmx_fio_fopen(fname, "w");
1119
1120     fprintf(out, "TITLE     %s\n", title);
1121     gmx_write_pdb_box(out, dd->bScrewPBC ? epbcSCREW : epbcXYZ, box);
1122     int molb = 0;
1123     for (int i = 0; i < natoms; i++)
1124     {
1125         int  ii = dd->globalAtomIndices[i];
1126         mtopGetAtomAndResidueName(mtop, ii, &molb, &atomname, &resnr, &resname, nullptr);
1127         int  c;
1128         real b;
1129         if (i < dd->comm->atomRanges.end(DDAtomRanges::Type::Zones))
1130         {
1131             c = 0;
1132             while (i >= dd->atomGrouping().subRange(0, dd->comm->zones.cg_range[c + 1]).end())
1133             {
1134                 c++;
1135             }
1136             b = c;
1137         }
1138         else if (i < dd->comm->atomRanges.end(DDAtomRanges::Type::Vsites))
1139         {
1140             b = dd->comm->zones.n;
1141         }
1142         else
1143         {
1144             b = dd->comm->zones.n + 1;
1145         }
1146         gmx_fprintf_pdb_atomline(out, epdbATOM, ii+1, atomname, ' ', resname, ' ', resnr, ' ',
1147                                  10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ], 1.0, b, "");
1148     }
1149     fprintf(out, "TER\n");
1150
1151     gmx_fio_fclose(out);
1152 }
1153
1154 real dd_cutoff_multibody(const gmx_domdec_t *dd)
1155 {
1156     gmx_domdec_comm_t *comm;
1157     int                di;
1158     real               r;
1159
1160     comm = dd->comm;
1161
1162     r = -1;
1163     if (comm->bInterCGBondeds)
1164     {
1165         if (comm->cutoff_mbody > 0)
1166         {
1167             r = comm->cutoff_mbody;
1168         }
1169         else
1170         {
1171             /* cutoff_mbody=0 means we do not have DLB */
1172             r = comm->cellsize_min[dd->dim[0]];
1173             for (di = 1; di < dd->ndim; di++)
1174             {
1175                 r = std::min(r, comm->cellsize_min[dd->dim[di]]);
1176             }
1177             if (comm->bBondComm)
1178             {
1179                 r = std::max(r, comm->cutoff_mbody);
1180             }
1181             else
1182             {
1183                 r = std::min(r, comm->cutoff);
1184             }
1185         }
1186     }
1187
1188     return r;
1189 }
1190
1191 real dd_cutoff_twobody(const gmx_domdec_t *dd)
1192 {
1193     real r_mb;
1194
1195     r_mb = dd_cutoff_multibody(dd);
1196
1197     return std::max(dd->comm->cutoff, r_mb);
1198 }
1199
1200
1201 static void dd_cart_coord2pmecoord(const gmx_domdec_t *dd, const ivec coord,
1202                                    ivec coord_pme)
1203 {
1204     int nc, ntot;
1205
1206     nc   = dd->nc[dd->comm->cartpmedim];
1207     ntot = dd->comm->ntot[dd->comm->cartpmedim];
1208     copy_ivec(coord, coord_pme);
1209     coord_pme[dd->comm->cartpmedim] =
1210         nc + (coord[dd->comm->cartpmedim]*(ntot - nc) + (ntot - nc)/2)/nc;
1211 }
1212
1213 static int ddindex2pmeindex(const gmx_domdec_t *dd, int ddindex)
1214 {
1215     int npp, npme;
1216
1217     npp  = dd->nnodes;
1218     npme = dd->comm->npmenodes;
1219
1220     /* Here we assign a PME node to communicate with this DD node
1221      * by assuming that the major index of both is x.
1222      * We add cr->npmenodes/2 to obtain an even distribution.
1223      */
1224     return (ddindex*npme + npme/2)/npp;
1225 }
1226
1227 static int *dd_interleaved_pme_ranks(const gmx_domdec_t *dd)
1228 {
1229     int *pme_rank;
1230     int  n, i, p0, p1;
1231
1232     snew(pme_rank, dd->comm->npmenodes);
1233     n = 0;
1234     for (i = 0; i < dd->nnodes; i++)
1235     {
1236         p0 = ddindex2pmeindex(dd, i);
1237         p1 = ddindex2pmeindex(dd, i+1);
1238         if (i+1 == dd->nnodes || p1 > p0)
1239         {
1240             if (debug)
1241             {
1242                 fprintf(debug, "pme_rank[%d] = %d\n", n, i+1+n);
1243             }
1244             pme_rank[n] = i + 1 + n;
1245             n++;
1246         }
1247     }
1248
1249     return pme_rank;
1250 }
1251
1252 static int gmx_ddcoord2pmeindex(const t_commrec *cr, int x, int y, int z)
1253 {
1254     gmx_domdec_t *dd;
1255     ivec          coords;
1256     int           slab;
1257
1258     dd = cr->dd;
1259     /*
1260        if (dd->comm->bCartesian) {
1261        gmx_ddindex2xyz(dd->nc,ddindex,coords);
1262        dd_coords2pmecoords(dd,coords,coords_pme);
1263        copy_ivec(dd->ntot,nc);
1264        nc[dd->cartpmedim]         -= dd->nc[dd->cartpmedim];
1265        coords_pme[dd->cartpmedim] -= dd->nc[dd->cartpmedim];
1266
1267        slab = (coords_pme[XX]*nc[YY] + coords_pme[YY])*nc[ZZ] + coords_pme[ZZ];
1268        } else {
1269        slab = (ddindex*cr->npmenodes + cr->npmenodes/2)/dd->nnodes;
1270        }
1271      */
1272     coords[XX] = x;
1273     coords[YY] = y;
1274     coords[ZZ] = z;
1275     slab       = ddindex2pmeindex(dd, dd_index(dd->nc, coords));
1276
1277     return slab;
1278 }
1279
1280 static int ddcoord2simnodeid(const t_commrec *cr, int x, int y, int z)
1281 {
1282     gmx_domdec_comm_t *comm;
1283     ivec               coords;
1284     int                ddindex, nodeid = -1;
1285
1286     comm = cr->dd->comm;
1287
1288     coords[XX] = x;
1289     coords[YY] = y;
1290     coords[ZZ] = z;
1291     if (comm->bCartesianPP_PME)
1292     {
1293 #if GMX_MPI
1294         MPI_Cart_rank(cr->mpi_comm_mysim, coords, &nodeid);
1295 #endif
1296     }
1297     else
1298     {
1299         ddindex = dd_index(cr->dd->nc, coords);
1300         if (comm->bCartesianPP)
1301         {
1302             nodeid = comm->ddindex2simnodeid[ddindex];
1303         }
1304         else
1305         {
1306             if (comm->pmenodes)
1307             {
1308                 nodeid = ddindex + gmx_ddcoord2pmeindex(cr, x, y, z);
1309             }
1310             else
1311             {
1312                 nodeid = ddindex;
1313             }
1314         }
1315     }
1316
1317     return nodeid;
1318 }
1319
1320 static int dd_simnode2pmenode(const gmx_domdec_t         *dd,
1321                               const t_commrec gmx_unused *cr,
1322                               int                         sim_nodeid)
1323 {
1324     int pmenode = -1;
1325
1326     const gmx_domdec_comm_t *comm = dd->comm;
1327
1328     /* This assumes a uniform x domain decomposition grid cell size */
1329     if (comm->bCartesianPP_PME)
1330     {
1331 #if GMX_MPI
1332         ivec coord, coord_pme;
1333         MPI_Cart_coords(cr->mpi_comm_mysim, sim_nodeid, DIM, coord);
1334         if (coord[comm->cartpmedim] < dd->nc[comm->cartpmedim])
1335         {
1336             /* This is a PP node */
1337             dd_cart_coord2pmecoord(dd, coord, coord_pme);
1338             MPI_Cart_rank(cr->mpi_comm_mysim, coord_pme, &pmenode);
1339         }
1340 #endif
1341     }
1342     else if (comm->bCartesianPP)
1343     {
1344         if (sim_nodeid < dd->nnodes)
1345         {
1346             pmenode = dd->nnodes + ddindex2pmeindex(dd, sim_nodeid);
1347         }
1348     }
1349     else
1350     {
1351         /* This assumes DD cells with identical x coordinates
1352          * are numbered sequentially.
1353          */
1354         if (dd->comm->pmenodes == nullptr)
1355         {
1356             if (sim_nodeid < dd->nnodes)
1357             {
1358                 /* The DD index equals the nodeid */
1359                 pmenode = dd->nnodes + ddindex2pmeindex(dd, sim_nodeid);
1360             }
1361         }
1362         else
1363         {
1364             int i = 0;
1365             while (sim_nodeid > dd->comm->pmenodes[i])
1366             {
1367                 i++;
1368             }
1369             if (sim_nodeid < dd->comm->pmenodes[i])
1370             {
1371                 pmenode = dd->comm->pmenodes[i];
1372             }
1373         }
1374     }
1375
1376     return pmenode;
1377 }
1378
1379 NumPmeDomains getNumPmeDomains(const gmx_domdec_t *dd)
1380 {
1381     if (dd != nullptr)
1382     {
1383         return {
1384                    dd->comm->npmenodes_x, dd->comm->npmenodes_y
1385         };
1386     }
1387     else
1388     {
1389         return {
1390                    1, 1
1391         };
1392     }
1393 }
1394
1395 std::vector<int> get_pme_ddranks(const t_commrec *cr, int pmenodeid)
1396 {
1397     gmx_domdec_t *dd;
1398     int           x, y, z;
1399     ivec          coord, coord_pme;
1400
1401     dd = cr->dd;
1402
1403     std::vector<int> ddranks;
1404     ddranks.reserve((dd->nnodes+cr->npmenodes-1)/cr->npmenodes);
1405
1406     for (x = 0; x < dd->nc[XX]; x++)
1407     {
1408         for (y = 0; y < dd->nc[YY]; y++)
1409         {
1410             for (z = 0; z < dd->nc[ZZ]; z++)
1411             {
1412                 if (dd->comm->bCartesianPP_PME)
1413                 {
1414                     coord[XX] = x;
1415                     coord[YY] = y;
1416                     coord[ZZ] = z;
1417                     dd_cart_coord2pmecoord(dd, coord, coord_pme);
1418                     if (dd->ci[XX] == coord_pme[XX] &&
1419                         dd->ci[YY] == coord_pme[YY] &&
1420                         dd->ci[ZZ] == coord_pme[ZZ])
1421                     {
1422                         ddranks.push_back(ddcoord2simnodeid(cr, x, y, z));
1423                     }
1424                 }
1425                 else
1426                 {
1427                     /* The slab corresponds to the nodeid in the PME group */
1428                     if (gmx_ddcoord2pmeindex(cr, x, y, z) == pmenodeid)
1429                     {
1430                         ddranks.push_back(ddcoord2simnodeid(cr, x, y, z));
1431                     }
1432                 }
1433             }
1434         }
1435     }
1436     return ddranks;
1437 }
1438
1439 static gmx_bool receive_vir_ener(const gmx_domdec_t *dd, const t_commrec *cr)
1440 {
1441     gmx_bool bReceive = TRUE;
1442
1443     if (cr->npmenodes < dd->nnodes)
1444     {
1445         gmx_domdec_comm_t *comm = dd->comm;
1446         if (comm->bCartesianPP_PME)
1447         {
1448 #if GMX_MPI
1449             int  pmenode = dd_simnode2pmenode(dd, cr, cr->sim_nodeid);
1450             ivec coords;
1451             MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, coords);
1452             coords[comm->cartpmedim]++;
1453             if (coords[comm->cartpmedim] < dd->nc[comm->cartpmedim])
1454             {
1455                 int rank;
1456                 MPI_Cart_rank(cr->mpi_comm_mysim, coords, &rank);
1457                 if (dd_simnode2pmenode(dd, cr, rank) == pmenode)
1458                 {
1459                     /* This is not the last PP node for pmenode */
1460                     bReceive = FALSE;
1461                 }
1462             }
1463 #else
1464             GMX_RELEASE_ASSERT(false, "Without MPI we should not have Cartesian PP-PME with #PMEnodes < #DDnodes");
1465 #endif
1466         }
1467         else
1468         {
1469             int pmenode = dd_simnode2pmenode(dd, cr, cr->sim_nodeid);
1470             if (cr->sim_nodeid+1 < cr->nnodes &&
1471                 dd_simnode2pmenode(dd, cr, cr->sim_nodeid+1) == pmenode)
1472             {
1473                 /* This is not the last PP node for pmenode */
1474                 bReceive = FALSE;
1475             }
1476         }
1477     }
1478
1479     return bReceive;
1480 }
1481
1482 static void set_zones_ncg_home(gmx_domdec_t *dd)
1483 {
1484     gmx_domdec_zones_t *zones;
1485     int                 i;
1486
1487     zones = &dd->comm->zones;
1488
1489     zones->cg_range[0] = 0;
1490     for (i = 1; i < zones->n+1; i++)
1491     {
1492         zones->cg_range[i] = dd->ncg_home;
1493     }
1494     /* zone_ncg1[0] should always be equal to ncg_home */
1495     dd->comm->zone_ncg1[0] = dd->ncg_home;
1496 }
1497
1498 static void restoreAtomGroups(gmx_domdec_t *dd,
1499                               const int *gcgs_index, const t_state *state)
1500 {
1501     gmx::ArrayRef<const int>  atomGroupsState        = state->cg_gl;
1502
1503     std::vector<int>         &globalAtomGroupIndices = dd->globalAtomGroupIndices;
1504     gmx::RangePartitioning   &atomGrouping           = dd->atomGrouping_;
1505
1506     globalAtomGroupIndices.resize(atomGroupsState.size());
1507     atomGrouping.clear();
1508
1509     /* Copy back the global charge group indices from state
1510      * and rebuild the local charge group to atom index.
1511      */
1512     for (gmx::index i = 0; i < atomGroupsState.size(); i++)
1513     {
1514         const int atomGroupGlobal  = atomGroupsState[i];
1515         const int groupSize        = gcgs_index[atomGroupGlobal + 1] - gcgs_index[atomGroupGlobal];
1516         globalAtomGroupIndices[i]  = atomGroupGlobal;
1517         atomGrouping.appendBlock(groupSize);
1518     }
1519
1520     dd->ncg_home = atomGroupsState.size();
1521     dd->comm->atomRanges.setEnd(DDAtomRanges::Type::Home, atomGrouping.fullRange().end());
1522
1523     set_zones_ncg_home(dd);
1524 }
1525
1526 static void dd_set_cginfo(gmx::ArrayRef<const int> index_gl, int cg0, int cg1,
1527                           t_forcerec *fr, char *bLocalCG)
1528 {
1529     cginfo_mb_t *cginfo_mb;
1530     int         *cginfo;
1531     int          cg;
1532
1533     if (fr != nullptr)
1534     {
1535         cginfo_mb = fr->cginfo_mb;
1536         cginfo    = fr->cginfo;
1537
1538         for (cg = cg0; cg < cg1; cg++)
1539         {
1540             cginfo[cg] = ddcginfo(cginfo_mb, index_gl[cg]);
1541         }
1542     }
1543
1544     if (bLocalCG != nullptr)
1545     {
1546         for (cg = cg0; cg < cg1; cg++)
1547         {
1548             bLocalCG[index_gl[cg]] = TRUE;
1549         }
1550     }
1551 }
1552
1553 static void make_dd_indices(gmx_domdec_t *dd,
1554                             const int *gcgs_index, int cg_start)
1555 {
1556     const int                 numZones               = dd->comm->zones.n;
1557     const int                *zone2cg                = dd->comm->zones.cg_range;
1558     const int                *zone_ncg1              = dd->comm->zone_ncg1;
1559     gmx::ArrayRef<const int>  globalAtomGroupIndices = dd->globalAtomGroupIndices;
1560     const gmx_bool            bCGs                   = dd->comm->bCGs;
1561
1562     std::vector<int>         &globalAtomIndices      = dd->globalAtomIndices;
1563
1564     if (zone2cg[1] != dd->ncg_home)
1565     {
1566         gmx_incons("dd->ncg_zone is not up to date");
1567     }
1568
1569     /* Make the local to global and global to local atom index */
1570     int a = dd->atomGrouping().subRange(cg_start, cg_start).begin();
1571     globalAtomIndices.resize(a);
1572     for (int zone = 0; zone < numZones; zone++)
1573     {
1574         int cg0;
1575         if (zone == 0)
1576         {
1577             cg0 = cg_start;
1578         }
1579         else
1580         {
1581             cg0 = zone2cg[zone];
1582         }
1583         int cg1    = zone2cg[zone+1];
1584         int cg1_p1 = cg0 + zone_ncg1[zone];
1585
1586         for (int cg = cg0; cg < cg1; cg++)
1587         {
1588             int zone1 = zone;
1589             if (cg >= cg1_p1)
1590             {
1591                 /* Signal that this cg is from more than one pulse away */
1592                 zone1 += numZones;
1593             }
1594             int cg_gl = globalAtomGroupIndices[cg];
1595             if (bCGs)
1596             {
1597                 for (int a_gl = gcgs_index[cg_gl]; a_gl < gcgs_index[cg_gl+1]; a_gl++)
1598                 {
1599                     globalAtomIndices.push_back(a_gl);
1600                     ga2la_set(dd->ga2la, a_gl, a, zone1);
1601                     a++;
1602                 }
1603             }
1604             else
1605             {
1606                 globalAtomIndices.push_back(cg_gl);
1607                 ga2la_set(dd->ga2la, cg_gl, a, zone1);
1608                 a++;
1609             }
1610         }
1611     }
1612 }
1613
1614 static int check_bLocalCG(gmx_domdec_t *dd, int ncg_sys, const char *bLocalCG,
1615                           const char *where)
1616 {
1617     int nerr = 0;
1618     if (bLocalCG == nullptr)
1619     {
1620         return nerr;
1621     }
1622     for (size_t i = 0; i < dd->globalAtomGroupIndices.size(); i++)
1623     {
1624         if (!bLocalCG[dd->globalAtomGroupIndices[i]])
1625         {
1626             fprintf(stderr,
1627                     "DD rank %d, %s: atom group %zu, global atom group %d is not marked in bLocalCG (ncg_home %d)\n", dd->rank, where, i + 1, dd->globalAtomGroupIndices[i] + 1, dd->ncg_home);
1628             nerr++;
1629         }
1630     }
1631     size_t ngl = 0;
1632     for (int i = 0; i < ncg_sys; i++)
1633     {
1634         if (bLocalCG[i])
1635         {
1636             ngl++;
1637         }
1638     }
1639     if (ngl != dd->globalAtomGroupIndices.size())
1640     {
1641         fprintf(stderr, "DD rank %d, %s: In bLocalCG %zu atom groups are marked as local, whereas there are %zu\n", dd->rank, where, ngl, dd->globalAtomGroupIndices.size());
1642         nerr++;
1643     }
1644
1645     return nerr;
1646 }
1647
1648 static void check_index_consistency(gmx_domdec_t *dd,
1649                                     int natoms_sys, int ncg_sys,
1650                                     const char *where)
1651 {
1652     int       nerr = 0;
1653
1654     const int numAtomsInZones = dd->comm->atomRanges.end(DDAtomRanges::Type::Zones);
1655
1656     if (dd->comm->DD_debug > 1)
1657     {
1658         std::vector<int> have(natoms_sys);
1659         for (int a = 0; a < numAtomsInZones; a++)
1660         {
1661             int globalAtomIndex = dd->globalAtomIndices[a];
1662             if (have[globalAtomIndex] > 0)
1663             {
1664                 fprintf(stderr, "DD rank %d: global atom %d occurs twice: index %d and %d\n", dd->rank, globalAtomIndex + 1, have[globalAtomIndex], a+1);
1665             }
1666             else
1667             {
1668                 have[globalAtomIndex] = a + 1;
1669             }
1670         }
1671     }
1672
1673     std::vector<int> have(numAtomsInZones);
1674
1675     int              ngl = 0;
1676     for (int i = 0; i < natoms_sys; i++)
1677     {
1678         int a;
1679         int cell;
1680         if (ga2la_get(dd->ga2la, i, &a, &cell))
1681         {
1682             if (a >= numAtomsInZones)
1683             {
1684                 fprintf(stderr, "DD rank %d: global atom %d marked as local atom %d, which is larger than nat_tot (%d)\n", dd->rank, i+1, a+1, numAtomsInZones);
1685                 nerr++;
1686             }
1687             else
1688             {
1689                 have[a] = 1;
1690                 if (dd->globalAtomIndices[a] != i)
1691                 {
1692                     fprintf(stderr, "DD rank %d: global atom %d marked as local atom %d, which has global atom index %d\n", dd->rank, i+1, a+1, dd->globalAtomIndices[a]+1);
1693                     nerr++;
1694                 }
1695             }
1696             ngl++;
1697         }
1698     }
1699     if (ngl != numAtomsInZones)
1700     {
1701         fprintf(stderr,
1702                 "DD rank %d, %s: %d global atom indices, %d local atoms\n",
1703                 dd->rank, where, ngl, numAtomsInZones);
1704     }
1705     for (int a = 0; a < numAtomsInZones; a++)
1706     {
1707         if (have[a] == 0)
1708         {
1709             fprintf(stderr,
1710                     "DD rank %d, %s: local atom %d, global %d has no global index\n",
1711                     dd->rank, where, a + 1, dd->globalAtomIndices[a] + 1);
1712         }
1713     }
1714
1715     nerr += check_bLocalCG(dd, ncg_sys, dd->comm->bLocalCG, where);
1716
1717     if (nerr > 0)
1718     {
1719         gmx_fatal(FARGS, "DD rank %d, %s: %d atom(group) index inconsistencies",
1720                   dd->rank, where, nerr);
1721     }
1722 }
1723
1724 /* Clear all DD global state indices, starting from \p atomGroupStart and \p atomStart */
1725 static void clearDDStateIndices(gmx_domdec_t *dd,
1726                                 int           atomGroupStart,
1727                                 int           atomStart)
1728 {
1729     if (atomStart == 0)
1730     {
1731         /* Clear the whole list without searching */
1732         ga2la_clear(dd->ga2la);
1733     }
1734     else
1735     {
1736         const int numAtomsInZones = dd->comm->atomRanges.end(DDAtomRanges::Type::Zones);
1737         for (int i = 0; i < numAtomsInZones; i++)
1738         {
1739             ga2la_del(dd->ga2la, dd->globalAtomIndices[i]);
1740         }
1741     }
1742
1743     char *bLocalCG = dd->comm->bLocalCG;
1744     if (bLocalCG)
1745     {
1746         for (size_t atomGroup = atomGroupStart; atomGroup < dd->globalAtomGroupIndices.size(); atomGroup++)
1747         {
1748             bLocalCG[dd->globalAtomGroupIndices[atomGroup]] = FALSE;
1749         }
1750     }
1751
1752     dd_clear_local_vsite_indices(dd);
1753
1754     if (dd->constraints)
1755     {
1756         dd_clear_local_constraint_indices(dd);
1757     }
1758 }
1759
1760 static bool check_grid_jump(int64_t             step,
1761                             const gmx_domdec_t *dd,
1762                             real                cutoff,
1763                             const gmx_ddbox_t  *ddbox,
1764                             gmx_bool            bFatal)
1765 {
1766     gmx_domdec_comm_t *comm    = dd->comm;
1767     bool               invalid = false;
1768
1769     for (int d = 1; d < dd->ndim; d++)
1770     {
1771         const DDCellsizesWithDlb &cellsizes = comm->cellsizesWithDlb[d];
1772         const int                 dim       = dd->dim[d];
1773         const real                limit     = grid_jump_limit(comm, cutoff, d);
1774         real                      bfac      = ddbox->box_size[dim];
1775         if (ddbox->tric_dir[dim])
1776         {
1777             bfac *= ddbox->skew_fac[dim];
1778         }
1779         if ((cellsizes.fracUpper - cellsizes.fracLowerMax)*bfac <  limit ||
1780             (cellsizes.fracLower - cellsizes.fracUpperMin)*bfac > -limit)
1781         {
1782             invalid = true;
1783
1784             if (bFatal)
1785             {
1786                 char buf[22];
1787
1788                 /* This error should never be triggered under normal
1789                  * circumstances, but you never know ...
1790                  */
1791                 gmx_fatal(FARGS, "step %s: The domain decomposition grid has shifted too much in the %c-direction around cell %d %d %d. This should not have happened. Running with fewer ranks might avoid this issue.",
1792                           gmx_step_str(step, buf),
1793                           dim2char(dim), dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
1794             }
1795         }
1796     }
1797
1798     return invalid;
1799 }
1800
1801 static float dd_force_load(gmx_domdec_comm_t *comm)
1802 {
1803     float load;
1804
1805     if (comm->eFlop)
1806     {
1807         load = comm->flop;
1808         if (comm->eFlop > 1)
1809         {
1810             load *= 1.0 + (comm->eFlop - 1)*(0.1*rand()/RAND_MAX - 0.05);
1811         }
1812     }
1813     else
1814     {
1815         load = comm->cycl[ddCyclF];
1816         if (comm->cycl_n[ddCyclF] > 1)
1817         {
1818             /* Subtract the maximum of the last n cycle counts
1819              * to get rid of possible high counts due to other sources,
1820              * for instance system activity, that would otherwise
1821              * affect the dynamic load balancing.
1822              */
1823             load -= comm->cycl_max[ddCyclF];
1824         }
1825
1826 #if GMX_MPI
1827         if (comm->cycl_n[ddCyclWaitGPU] && comm->nrank_gpu_shared > 1)
1828         {
1829             float gpu_wait, gpu_wait_sum;
1830
1831             gpu_wait = comm->cycl[ddCyclWaitGPU];
1832             if (comm->cycl_n[ddCyclF] > 1)
1833             {
1834                 /* We should remove the WaitGPU time of the same MD step
1835                  * as the one with the maximum F time, since the F time
1836                  * and the wait time are not independent.
1837                  * Furthermore, the step for the max F time should be chosen
1838                  * the same on all ranks that share the same GPU.
1839                  * But to keep the code simple, we remove the average instead.
1840                  * The main reason for artificially long times at some steps
1841                  * is spurious CPU activity or MPI time, so we don't expect
1842                  * that changes in the GPU wait time matter a lot here.
1843                  */
1844                 gpu_wait *= (comm->cycl_n[ddCyclF] - 1)/(float)comm->cycl_n[ddCyclF];
1845             }
1846             /* Sum the wait times over the ranks that share the same GPU */
1847             MPI_Allreduce(&gpu_wait, &gpu_wait_sum, 1, MPI_FLOAT, MPI_SUM,
1848                           comm->mpi_comm_gpu_shared);
1849             /* Replace the wait time by the average over the ranks */
1850             load += -gpu_wait + gpu_wait_sum/comm->nrank_gpu_shared;
1851         }
1852 #endif
1853     }
1854
1855     return load;
1856 }
1857
1858 static void set_slb_pme_dim_f(gmx_domdec_t *dd, int dim, real **dim_f)
1859 {
1860     gmx_domdec_comm_t *comm;
1861     int                i;
1862
1863     comm = dd->comm;
1864
1865     snew(*dim_f, dd->nc[dim]+1);
1866     (*dim_f)[0] = 0;
1867     for (i = 1; i < dd->nc[dim]; i++)
1868     {
1869         if (comm->slb_frac[dim])
1870         {
1871             (*dim_f)[i] = (*dim_f)[i-1] + comm->slb_frac[dim][i-1];
1872         }
1873         else
1874         {
1875             (*dim_f)[i] = (real)i/(real)dd->nc[dim];
1876         }
1877     }
1878     (*dim_f)[dd->nc[dim]] = 1;
1879 }
1880
1881 static void init_ddpme(gmx_domdec_t *dd, gmx_ddpme_t *ddpme, int dimind)
1882 {
1883     int  pmeindex, slab, nso, i;
1884     ivec xyz;
1885
1886     if (dimind == 0 && dd->dim[0] == YY && dd->comm->npmenodes_x == 1)
1887     {
1888         ddpme->dim = YY;
1889     }
1890     else
1891     {
1892         ddpme->dim = dimind;
1893     }
1894     ddpme->dim_match = (ddpme->dim == dd->dim[dimind]);
1895
1896     ddpme->nslab = (ddpme->dim == 0 ?
1897                     dd->comm->npmenodes_x :
1898                     dd->comm->npmenodes_y);
1899
1900     if (ddpme->nslab <= 1)
1901     {
1902         return;
1903     }
1904
1905     nso = dd->comm->npmenodes/ddpme->nslab;
1906     /* Determine for each PME slab the PP location range for dimension dim */
1907     snew(ddpme->pp_min, ddpme->nslab);
1908     snew(ddpme->pp_max, ddpme->nslab);
1909     for (slab = 0; slab < ddpme->nslab; slab++)
1910     {
1911         ddpme->pp_min[slab] = dd->nc[dd->dim[dimind]] - 1;
1912         ddpme->pp_max[slab] = 0;
1913     }
1914     for (i = 0; i < dd->nnodes; i++)
1915     {
1916         ddindex2xyz(dd->nc, i, xyz);
1917         /* For y only use our y/z slab.
1918          * This assumes that the PME x grid size matches the DD grid size.
1919          */
1920         if (dimind == 0 || xyz[XX] == dd->ci[XX])
1921         {
1922             pmeindex = ddindex2pmeindex(dd, i);
1923             if (dimind == 0)
1924             {
1925                 slab = pmeindex/nso;
1926             }
1927             else
1928             {
1929                 slab = pmeindex % ddpme->nslab;
1930             }
1931             ddpme->pp_min[slab] = std::min(ddpme->pp_min[slab], xyz[dimind]);
1932             ddpme->pp_max[slab] = std::max(ddpme->pp_max[slab], xyz[dimind]);
1933         }
1934     }
1935
1936     set_slb_pme_dim_f(dd, ddpme->dim, &ddpme->slb_dim_f);
1937 }
1938
1939 int dd_pme_maxshift_x(const gmx_domdec_t *dd)
1940 {
1941     if (dd->comm->ddpme[0].dim == XX)
1942     {
1943         return dd->comm->ddpme[0].maxshift;
1944     }
1945     else
1946     {
1947         return 0;
1948     }
1949 }
1950
1951 int dd_pme_maxshift_y(const gmx_domdec_t *dd)
1952 {
1953     if (dd->comm->ddpme[0].dim == YY)
1954     {
1955         return dd->comm->ddpme[0].maxshift;
1956     }
1957     else if (dd->comm->npmedecompdim >= 2 && dd->comm->ddpme[1].dim == YY)
1958     {
1959         return dd->comm->ddpme[1].maxshift;
1960     }
1961     else
1962     {
1963         return 0;
1964     }
1965 }
1966
1967 static void comm_dd_ns_cell_sizes(gmx_domdec_t *dd,
1968                                   gmx_ddbox_t *ddbox,
1969                                   rvec cell_ns_x0, rvec cell_ns_x1,
1970                                   int64_t step)
1971 {
1972     gmx_domdec_comm_t *comm;
1973     int                dim_ind, dim;
1974
1975     comm = dd->comm;
1976
1977     for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
1978     {
1979         dim = dd->dim[dim_ind];
1980
1981         /* Without PBC we don't have restrictions on the outer cells */
1982         if (!(dim >= ddbox->npbcdim &&
1983               (dd->ci[dim] == 0 || dd->ci[dim] == dd->nc[dim] - 1)) &&
1984             isDlbOn(comm) &&
1985             (comm->cell_x1[dim] - comm->cell_x0[dim])*ddbox->skew_fac[dim] <
1986             comm->cellsize_min[dim])
1987         {
1988             char buf[22];
1989             gmx_fatal(FARGS, "step %s: The %c-size (%f) times the triclinic skew factor (%f) is smaller than the smallest allowed cell size (%f) for domain decomposition grid cell %d %d %d",
1990                       gmx_step_str(step, buf), dim2char(dim),
1991                       comm->cell_x1[dim] - comm->cell_x0[dim],
1992                       ddbox->skew_fac[dim],
1993                       dd->comm->cellsize_min[dim],
1994                       dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
1995         }
1996     }
1997
1998     if ((isDlbOn(dd->comm) && dd->ndim > 1) || ddbox->nboundeddim < DIM)
1999     {
2000         /* Communicate the boundaries and update cell_ns_x0/1 */
2001         dd_move_cellx(dd, ddbox, cell_ns_x0, cell_ns_x1);
2002         if (isDlbOn(dd->comm) && dd->ndim > 1)
2003         {
2004             check_grid_jump(step, dd, dd->comm->cutoff, ddbox, TRUE);
2005         }
2006     }
2007 }
2008
2009 void dd_cycles_add(const gmx_domdec_t *dd, float cycles, int ddCycl)
2010 {
2011     /* Note that the cycles value can be incorrect, either 0 or some
2012      * extremely large value, when our thread migrated to another core
2013      * with an unsynchronized cycle counter. If this happens less often
2014      * that once per nstlist steps, this will not cause issues, since
2015      * we later subtract the maximum value from the sum over nstlist steps.
2016      * A zero count will slightly lower the total, but that's a small effect.
2017      * Note that the main purpose of the subtraction of the maximum value
2018      * is to avoid throwing off the load balancing when stalls occur due
2019      * e.g. system activity or network congestion.
2020      */
2021     dd->comm->cycl[ddCycl] += cycles;
2022     dd->comm->cycl_n[ddCycl]++;
2023     if (cycles > dd->comm->cycl_max[ddCycl])
2024     {
2025         dd->comm->cycl_max[ddCycl] = cycles;
2026     }
2027 }
2028
2029 static double force_flop_count(t_nrnb *nrnb)
2030 {
2031     int         i;
2032     double      sum;
2033     const char *name;
2034
2035     sum = 0;
2036     for (i = 0; i < eNR_NBKERNEL_FREE_ENERGY; i++)
2037     {
2038         /* To get closer to the real timings, we half the count
2039          * for the normal loops and again half it for water loops.
2040          */
2041         name = nrnb_str(i);
2042         if (strstr(name, "W3") != nullptr || strstr(name, "W4") != nullptr)
2043         {
2044             sum += nrnb->n[i]*0.25*cost_nrnb(i);
2045         }
2046         else
2047         {
2048             sum += nrnb->n[i]*0.50*cost_nrnb(i);
2049         }
2050     }
2051     for (i = eNR_NBKERNEL_FREE_ENERGY; i <= eNR_NB14; i++)
2052     {
2053         name = nrnb_str(i);
2054         if (strstr(name, "W3") != nullptr || strstr(name, "W4") != nullptr)
2055         {
2056             sum += nrnb->n[i]*cost_nrnb(i);
2057         }
2058     }
2059     for (i = eNR_BONDS; i <= eNR_WALLS; i++)
2060     {
2061         sum += nrnb->n[i]*cost_nrnb(i);
2062     }
2063
2064     return sum;
2065 }
2066
2067 void dd_force_flop_start(gmx_domdec_t *dd, t_nrnb *nrnb)
2068 {
2069     if (dd->comm->eFlop)
2070     {
2071         dd->comm->flop -= force_flop_count(nrnb);
2072     }
2073 }
2074 void dd_force_flop_stop(gmx_domdec_t *dd, t_nrnb *nrnb)
2075 {
2076     if (dd->comm->eFlop)
2077     {
2078         dd->comm->flop += force_flop_count(nrnb);
2079         dd->comm->flop_n++;
2080     }
2081 }
2082
2083 static void clear_dd_cycle_counts(gmx_domdec_t *dd)
2084 {
2085     int i;
2086
2087     for (i = 0; i < ddCyclNr; i++)
2088     {
2089         dd->comm->cycl[i]     = 0;
2090         dd->comm->cycl_n[i]   = 0;
2091         dd->comm->cycl_max[i] = 0;
2092     }
2093     dd->comm->flop   = 0;
2094     dd->comm->flop_n = 0;
2095 }
2096
2097 static void get_load_distribution(gmx_domdec_t *dd, gmx_wallcycle_t wcycle)
2098 {
2099     gmx_domdec_comm_t *comm;
2100     domdec_load_t     *load;
2101     float              cell_frac = 0, sbuf[DD_NLOAD_MAX];
2102     gmx_bool           bSepPME;
2103
2104     if (debug)
2105     {
2106         fprintf(debug, "get_load_distribution start\n");
2107     }
2108
2109     wallcycle_start(wcycle, ewcDDCOMMLOAD);
2110
2111     comm = dd->comm;
2112
2113     bSepPME = (dd->pme_nodeid >= 0);
2114
2115     if (dd->ndim == 0 && bSepPME)
2116     {
2117         /* Without decomposition, but with PME nodes, we need the load */
2118         comm->load[0].mdf = comm->cycl[ddCyclPPduringPME];
2119         comm->load[0].pme = comm->cycl[ddCyclPME];
2120     }
2121
2122     for (int d = dd->ndim - 1; d >= 0; d--)
2123     {
2124         const DDCellsizesWithDlb &cellsizes = comm->cellsizesWithDlb[d];
2125         const int                 dim       = dd->dim[d];
2126         /* Check if we participate in the communication in this dimension */
2127         if (d == dd->ndim-1 ||
2128             (dd->ci[dd->dim[d+1]] == 0 && dd->ci[dd->dim[dd->ndim-1]] == 0))
2129         {
2130             load = &comm->load[d];
2131             if (isDlbOn(dd->comm))
2132             {
2133                 cell_frac = cellsizes.fracUpper - cellsizes.fracLower;
2134             }
2135             int pos = 0;
2136             if (d == dd->ndim-1)
2137             {
2138                 sbuf[pos++] = dd_force_load(comm);
2139                 sbuf[pos++] = sbuf[0];
2140                 if (isDlbOn(dd->comm))
2141                 {
2142                     sbuf[pos++] = sbuf[0];
2143                     sbuf[pos++] = cell_frac;
2144                     if (d > 0)
2145                     {
2146                         sbuf[pos++] = cellsizes.fracLowerMax;
2147                         sbuf[pos++] = cellsizes.fracUpperMin;
2148                     }
2149                 }
2150                 if (bSepPME)
2151                 {
2152                     sbuf[pos++] = comm->cycl[ddCyclPPduringPME];
2153                     sbuf[pos++] = comm->cycl[ddCyclPME];
2154                 }
2155             }
2156             else
2157             {
2158                 sbuf[pos++] = comm->load[d+1].sum;
2159                 sbuf[pos++] = comm->load[d+1].max;
2160                 if (isDlbOn(dd->comm))
2161                 {
2162                     sbuf[pos++] = comm->load[d+1].sum_m;
2163                     sbuf[pos++] = comm->load[d+1].cvol_min*cell_frac;
2164                     sbuf[pos++] = comm->load[d+1].flags;
2165                     if (d > 0)
2166                     {
2167                         sbuf[pos++] = cellsizes.fracLowerMax;
2168                         sbuf[pos++] = cellsizes.fracUpperMin;
2169                     }
2170                 }
2171                 if (bSepPME)
2172                 {
2173                     sbuf[pos++] = comm->load[d+1].mdf;
2174                     sbuf[pos++] = comm->load[d+1].pme;
2175                 }
2176             }
2177             load->nload = pos;
2178             /* Communicate a row in DD direction d.
2179              * The communicators are setup such that the root always has rank 0.
2180              */
2181 #if GMX_MPI
2182             MPI_Gather(sbuf, load->nload*sizeof(float), MPI_BYTE,
2183                        load->load, load->nload*sizeof(float), MPI_BYTE,
2184                        0, comm->mpi_comm_load[d]);
2185 #endif
2186             if (dd->ci[dim] == dd->master_ci[dim])
2187             {
2188                 /* We are the master along this row, process this row */
2189                 RowMaster *rowMaster = nullptr;
2190
2191                 if (isDlbOn(comm))
2192                 {
2193                     rowMaster = cellsizes.rowMaster.get();
2194                 }
2195                 load->sum      = 0;
2196                 load->max      = 0;
2197                 load->sum_m    = 0;
2198                 load->cvol_min = 1;
2199                 load->flags    = 0;
2200                 load->mdf      = 0;
2201                 load->pme      = 0;
2202                 int pos        = 0;
2203                 for (int i = 0; i < dd->nc[dim]; i++)
2204                 {
2205                     load->sum += load->load[pos++];
2206                     load->max  = std::max(load->max, load->load[pos]);
2207                     pos++;
2208                     if (isDlbOn(dd->comm))
2209                     {
2210                         if (rowMaster->dlbIsLimited)
2211                         {
2212                             /* This direction could not be load balanced properly,
2213                              * therefore we need to use the maximum iso the average load.
2214                              */
2215                             load->sum_m = std::max(load->sum_m, load->load[pos]);
2216                         }
2217                         else
2218                         {
2219                             load->sum_m += load->load[pos];
2220                         }
2221                         pos++;
2222                         load->cvol_min = std::min(load->cvol_min, load->load[pos]);
2223                         pos++;
2224                         if (d < dd->ndim-1)
2225                         {
2226                             load->flags = (int)(load->load[pos++] + 0.5);
2227                         }
2228                         if (d > 0)
2229                         {
2230                             rowMaster->bounds[i].cellFracLowerMax = load->load[pos++];
2231                             rowMaster->bounds[i].cellFracUpperMin = load->load[pos++];
2232                         }
2233                     }
2234                     if (bSepPME)
2235                     {
2236                         load->mdf = std::max(load->mdf, load->load[pos]);
2237                         pos++;
2238                         load->pme = std::max(load->pme, load->load[pos]);
2239                         pos++;
2240                     }
2241                 }
2242                 if (isDlbOn(comm) && rowMaster->dlbIsLimited)
2243                 {
2244                     load->sum_m *= dd->nc[dim];
2245                     load->flags |= (1<<d);
2246                 }
2247             }
2248         }
2249     }
2250
2251     if (DDMASTER(dd))
2252     {
2253         comm->nload      += dd_load_count(comm);
2254         comm->load_step  += comm->cycl[ddCyclStep];
2255         comm->load_sum   += comm->load[0].sum;
2256         comm->load_max   += comm->load[0].max;
2257         if (isDlbOn(comm))
2258         {
2259             for (int d = 0; d < dd->ndim; d++)
2260             {
2261                 if (comm->load[0].flags & (1<<d))
2262                 {
2263                     comm->load_lim[d]++;
2264                 }
2265             }
2266         }
2267         if (bSepPME)
2268         {
2269             comm->load_mdf += comm->load[0].mdf;
2270             comm->load_pme += comm->load[0].pme;
2271         }
2272     }
2273
2274     wallcycle_stop(wcycle, ewcDDCOMMLOAD);
2275
2276     if (debug)
2277     {
2278         fprintf(debug, "get_load_distribution finished\n");
2279     }
2280 }
2281
2282 static float dd_force_load_fraction(gmx_domdec_t *dd)
2283 {
2284     /* Return the relative performance loss on the total run time
2285      * due to the force calculation load imbalance.
2286      */
2287     if (dd->comm->nload > 0 && dd->comm->load_step > 0)
2288     {
2289         return dd->comm->load_sum/(dd->comm->load_step*dd->nnodes);
2290     }
2291     else
2292     {
2293         return 0;
2294     }
2295 }
2296
2297 static float dd_force_imb_perf_loss(gmx_domdec_t *dd)
2298 {
2299     /* Return the relative performance loss on the total run time
2300      * due to the force calculation load imbalance.
2301      */
2302     if (dd->comm->nload > 0 && dd->comm->load_step > 0)
2303     {
2304         return
2305             (dd->comm->load_max*dd->nnodes - dd->comm->load_sum)/
2306             (dd->comm->load_step*dd->nnodes);
2307     }
2308     else
2309     {
2310         return 0;
2311     }
2312 }
2313
2314 static void print_dd_load_av(FILE *fplog, gmx_domdec_t *dd)
2315 {
2316     gmx_domdec_comm_t *comm = dd->comm;
2317
2318     /* Only the master rank prints loads and only if we measured loads */
2319     if (!DDMASTER(dd) || comm->nload == 0)
2320     {
2321         return;
2322     }
2323
2324     char  buf[STRLEN];
2325     int   numPpRanks   = dd->nnodes;
2326     int   numPmeRanks  = (dd->pme_nodeid >= 0) ? comm->npmenodes : 0;
2327     int   numRanks     = numPpRanks + numPmeRanks;
2328     float lossFraction = 0;
2329
2330     /* Print the average load imbalance and performance loss */
2331     if (dd->nnodes > 1 && comm->load_sum > 0)
2332     {
2333         float imbalance = comm->load_max*numPpRanks/comm->load_sum - 1;
2334         lossFraction    = dd_force_imb_perf_loss(dd);
2335
2336         std::string msg         = "\n Dynamic load balancing report:\n";
2337         std::string dlbStateStr;
2338
2339         switch (dd->comm->dlbState)
2340         {
2341             case edlbsOffUser:
2342                 dlbStateStr = "DLB was off during the run per user request.";
2343                 break;
2344             case edlbsOffForever:
2345                 /* Currectly this can happen due to performance loss observed, cell size
2346                  * limitations or incompatibility with other settings observed during
2347                  * determineInitialDlbState(). */
2348                 dlbStateStr = "DLB got disabled because it was unsuitable to use.";
2349                 break;
2350             case edlbsOffCanTurnOn:
2351                 dlbStateStr = "DLB was off during the run due to low measured imbalance.";
2352                 break;
2353             case edlbsOffTemporarilyLocked:
2354                 dlbStateStr = "DLB was locked at the end of the run due to unfinished PP-PME balancing.";
2355                 break;
2356             case edlbsOnCanTurnOff:
2357                 dlbStateStr = "DLB was turned on during the run due to measured imbalance.";
2358                 break;
2359             case edlbsOnUser:
2360                 dlbStateStr = "DLB was permanently on during the run per user request.";
2361                 break;
2362             default:
2363                 GMX_ASSERT(false, "Undocumented DLB state");
2364         }
2365
2366         msg += " " + dlbStateStr + "\n";
2367         msg += gmx::formatString(" Average load imbalance: %.1f%%.\n", imbalance*100);
2368         msg += gmx::formatString(" The balanceable part of the MD step is %d%%, load imbalance is computed from this.\n",
2369                                  static_cast<int>(dd_force_load_fraction(dd)*100 + 0.5));
2370         msg += gmx::formatString(" Part of the total run time spent waiting due to load imbalance: %.1f%%.\n",
2371                                  lossFraction*100);
2372         fprintf(fplog, "%s", msg.c_str());
2373         fprintf(stderr, "%s", msg.c_str());
2374     }
2375
2376     /* Print during what percentage of steps the  load balancing was limited */
2377     bool dlbWasLimited = false;
2378     if (isDlbOn(comm))
2379     {
2380         sprintf(buf, " Steps where the load balancing was limited by -rdd, -rcon and/or -dds:");
2381         for (int d = 0; d < dd->ndim; d++)
2382         {
2383             int limitPercentage = (200*comm->load_lim[d] + 1)/(2*comm->nload);
2384             sprintf(buf+strlen(buf), " %c %d %%",
2385                     dim2char(dd->dim[d]), limitPercentage);
2386             if (limitPercentage >= 50)
2387             {
2388                 dlbWasLimited = true;
2389             }
2390         }
2391         sprintf(buf + strlen(buf), "\n");
2392         fprintf(fplog, "%s", buf);
2393         fprintf(stderr, "%s", buf);
2394     }
2395
2396     /* Print the performance loss due to separate PME - PP rank imbalance */
2397     float lossFractionPme = 0;
2398     if (numPmeRanks > 0 && comm->load_mdf > 0 && comm->load_step > 0)
2399     {
2400         float pmeForceRatio = comm->load_pme/comm->load_mdf;
2401         lossFractionPme     = (comm->load_pme - comm->load_mdf)/comm->load_step;
2402         if (lossFractionPme <= 0)
2403         {
2404             lossFractionPme *= numPmeRanks/static_cast<float>(numRanks);
2405         }
2406         else
2407         {
2408             lossFractionPme *= numPpRanks/static_cast<float>(numRanks);
2409         }
2410         sprintf(buf, " Average PME mesh/force load: %5.3f\n", pmeForceRatio);
2411         fprintf(fplog, "%s", buf);
2412         fprintf(stderr, "%s", buf);
2413         sprintf(buf, " Part of the total run time spent waiting due to PP/PME imbalance: %.1f %%\n", std::fabs(lossFractionPme)*100);
2414         fprintf(fplog, "%s", buf);
2415         fprintf(stderr, "%s", buf);
2416     }
2417     fprintf(fplog, "\n");
2418     fprintf(stderr, "\n");
2419
2420     if (lossFraction >= DD_PERF_LOSS_WARN)
2421     {
2422         sprintf(buf,
2423                 "NOTE: %.1f %% of the available CPU time was lost due to load imbalance\n"
2424                 "      in the domain decomposition.\n", lossFraction*100);
2425         if (!isDlbOn(comm))
2426         {
2427             sprintf(buf+strlen(buf), "      You might want to use dynamic load balancing (option -dlb.)\n");
2428         }
2429         else if (dlbWasLimited)
2430         {
2431             sprintf(buf+strlen(buf), "      You might want to decrease the cell size limit (options -rdd, -rcon and/or -dds).\n");
2432         }
2433         fprintf(fplog, "%s\n", buf);
2434         fprintf(stderr, "%s\n", buf);
2435     }
2436     if (numPmeRanks > 0 && std::fabs(lossFractionPme) >= DD_PERF_LOSS_WARN)
2437     {
2438         sprintf(buf,
2439                 "NOTE: %.1f %% performance was lost because the PME ranks\n"
2440                 "      had %s work to do than the PP ranks.\n"
2441                 "      You might want to %s the number of PME ranks\n"
2442                 "      or %s the cut-off and the grid spacing.\n",
2443                 std::fabs(lossFractionPme*100),
2444                 (lossFractionPme < 0) ? "less"     : "more",
2445                 (lossFractionPme < 0) ? "decrease" : "increase",
2446                 (lossFractionPme < 0) ? "decrease" : "increase");
2447         fprintf(fplog, "%s\n", buf);
2448         fprintf(stderr, "%s\n", buf);
2449     }
2450 }
2451
2452 static float dd_vol_min(gmx_domdec_t *dd)
2453 {
2454     return dd->comm->load[0].cvol_min*dd->nnodes;
2455 }
2456
2457 static gmx_bool dd_load_flags(gmx_domdec_t *dd)
2458 {
2459     return dd->comm->load[0].flags;
2460 }
2461
2462 static float dd_f_imbal(gmx_domdec_t *dd)
2463 {
2464     if (dd->comm->load[0].sum > 0)
2465     {
2466         return dd->comm->load[0].max*dd->nnodes/dd->comm->load[0].sum - 1.0f;
2467     }
2468     else
2469     {
2470         /* Something is wrong in the cycle counting, report no load imbalance */
2471         return 0.0f;
2472     }
2473 }
2474
2475 float dd_pme_f_ratio(gmx_domdec_t *dd)
2476 {
2477     /* Should only be called on the DD master rank */
2478     assert(DDMASTER(dd));
2479
2480     if (dd->comm->load[0].mdf > 0 && dd->comm->cycl_n[ddCyclPME] > 0)
2481     {
2482         return dd->comm->load[0].pme/dd->comm->load[0].mdf;
2483     }
2484     else
2485     {
2486         return -1.0;
2487     }
2488 }
2489
2490 static void dd_print_load(FILE *fplog, gmx_domdec_t *dd, int64_t step)
2491 {
2492     int  flags, d;
2493     char buf[22];
2494
2495     flags = dd_load_flags(dd);
2496     if (flags)
2497     {
2498         fprintf(fplog,
2499                 "DD  load balancing is limited by minimum cell size in dimension");
2500         for (d = 0; d < dd->ndim; d++)
2501         {
2502             if (flags & (1<<d))
2503             {
2504                 fprintf(fplog, " %c", dim2char(dd->dim[d]));
2505             }
2506         }
2507         fprintf(fplog, "\n");
2508     }
2509     fprintf(fplog, "DD  step %s", gmx_step_str(step, buf));
2510     if (isDlbOn(dd->comm))
2511     {
2512         fprintf(fplog, "  vol min/aver %5.3f%c",
2513                 dd_vol_min(dd), flags ? '!' : ' ');
2514     }
2515     if (dd->nnodes > 1)
2516     {
2517         fprintf(fplog, " load imb.: force %4.1f%%", dd_f_imbal(dd)*100);
2518     }
2519     if (dd->comm->cycl_n[ddCyclPME])
2520     {
2521         fprintf(fplog, "  pme mesh/force %5.3f", dd_pme_f_ratio(dd));
2522     }
2523     fprintf(fplog, "\n\n");
2524 }
2525
2526 static void dd_print_load_verbose(gmx_domdec_t *dd)
2527 {
2528     if (isDlbOn(dd->comm))
2529     {
2530         fprintf(stderr, "vol %4.2f%c ",
2531                 dd_vol_min(dd), dd_load_flags(dd) ? '!' : ' ');
2532     }
2533     if (dd->nnodes > 1)
2534     {
2535         fprintf(stderr, "imb F %2d%% ", (int)(dd_f_imbal(dd)*100+0.5));
2536     }
2537     if (dd->comm->cycl_n[ddCyclPME])
2538     {
2539         fprintf(stderr, "pme/F %4.2f ", dd_pme_f_ratio(dd));
2540     }
2541 }
2542
2543 #if GMX_MPI
2544 static void make_load_communicator(gmx_domdec_t *dd, int dim_ind, ivec loc)
2545 {
2546     MPI_Comm           c_row;
2547     int                dim, i, rank;
2548     ivec               loc_c;
2549     gmx_bool           bPartOfGroup = FALSE;
2550
2551     dim = dd->dim[dim_ind];
2552     copy_ivec(loc, loc_c);
2553     for (i = 0; i < dd->nc[dim]; i++)
2554     {
2555         loc_c[dim] = i;
2556         rank       = dd_index(dd->nc, loc_c);
2557         if (rank == dd->rank)
2558         {
2559             /* This process is part of the group */
2560             bPartOfGroup = TRUE;
2561         }
2562     }
2563     MPI_Comm_split(dd->mpi_comm_all, bPartOfGroup ? 0 : MPI_UNDEFINED, dd->rank,
2564                    &c_row);
2565     if (bPartOfGroup)
2566     {
2567         dd->comm->mpi_comm_load[dim_ind] = c_row;
2568         if (!isDlbDisabled(dd->comm))
2569         {
2570             DDCellsizesWithDlb &cellsizes = dd->comm->cellsizesWithDlb[dim_ind];
2571
2572             if (dd->ci[dim] == dd->master_ci[dim])
2573             {
2574                 /* This is the root process of this row */
2575                 cellsizes.rowMaster  = gmx::compat::make_unique<RowMaster>();
2576
2577                 RowMaster &rowMaster = *cellsizes.rowMaster;
2578                 rowMaster.cellFrac.resize(ddCellFractionBufferSize(dd, dim_ind));
2579                 rowMaster.oldCellFrac.resize(dd->nc[dim] + 1);
2580                 rowMaster.isCellMin.resize(dd->nc[dim]);
2581                 if (dim_ind > 0)
2582                 {
2583                     rowMaster.bounds.resize(dd->nc[dim]);
2584                 }
2585                 rowMaster.buf_ncd.resize(dd->nc[dim]);
2586             }
2587             else
2588             {
2589                 /* This is not a root process, we only need to receive cell_f */
2590                 cellsizes.fracRow.resize(ddCellFractionBufferSize(dd, dim_ind));
2591             }
2592         }
2593         if (dd->ci[dim] == dd->master_ci[dim])
2594         {
2595             snew(dd->comm->load[dim_ind].load, dd->nc[dim]*DD_NLOAD_MAX);
2596         }
2597     }
2598 }
2599 #endif
2600
2601 void dd_setup_dlb_resource_sharing(t_commrec            *cr,
2602                                    int                   gpu_id)
2603 {
2604 #if GMX_MPI
2605     int           physicalnode_id_hash;
2606     gmx_domdec_t *dd;
2607     MPI_Comm      mpi_comm_pp_physicalnode;
2608
2609     if (!thisRankHasDuty(cr, DUTY_PP) || gpu_id < 0)
2610     {
2611         /* Only ranks with short-ranged tasks (currently) use GPUs.
2612          * If we don't have GPUs assigned, there are no resources to share.
2613          */
2614         return;
2615     }
2616
2617     physicalnode_id_hash = gmx_physicalnode_id_hash();
2618
2619     dd = cr->dd;
2620
2621     if (debug)
2622     {
2623         fprintf(debug, "dd_setup_dd_dlb_gpu_sharing:\n");
2624         fprintf(debug, "DD PP rank %d physical node hash %d gpu_id %d\n",
2625                 dd->rank, physicalnode_id_hash, gpu_id);
2626     }
2627     /* Split the PP communicator over the physical nodes */
2628     /* TODO: See if we should store this (before), as it's also used for
2629      * for the nodecomm summation.
2630      */
2631     // TODO PhysicalNodeCommunicator could be extended/used to handle
2632     // the need for per-node per-group communicators.
2633     MPI_Comm_split(dd->mpi_comm_all, physicalnode_id_hash, dd->rank,
2634                    &mpi_comm_pp_physicalnode);
2635     MPI_Comm_split(mpi_comm_pp_physicalnode, gpu_id, dd->rank,
2636                    &dd->comm->mpi_comm_gpu_shared);
2637     MPI_Comm_free(&mpi_comm_pp_physicalnode);
2638     MPI_Comm_size(dd->comm->mpi_comm_gpu_shared, &dd->comm->nrank_gpu_shared);
2639
2640     if (debug)
2641     {
2642         fprintf(debug, "nrank_gpu_shared %d\n", dd->comm->nrank_gpu_shared);
2643     }
2644
2645     /* Note that some ranks could share a GPU, while others don't */
2646
2647     if (dd->comm->nrank_gpu_shared == 1)
2648     {
2649         MPI_Comm_free(&dd->comm->mpi_comm_gpu_shared);
2650     }
2651 #else
2652     GMX_UNUSED_VALUE(cr);
2653     GMX_UNUSED_VALUE(gpu_id);
2654 #endif
2655 }
2656
2657 static void make_load_communicators(gmx_domdec_t gmx_unused *dd)
2658 {
2659 #if GMX_MPI
2660     int  dim0, dim1, i, j;
2661     ivec loc;
2662
2663     if (debug)
2664     {
2665         fprintf(debug, "Making load communicators\n");
2666     }
2667
2668     snew(dd->comm->load,          std::max(dd->ndim, 1));
2669     snew(dd->comm->mpi_comm_load, std::max(dd->ndim, 1));
2670
2671     if (dd->ndim == 0)
2672     {
2673         return;
2674     }
2675
2676     clear_ivec(loc);
2677     make_load_communicator(dd, 0, loc);
2678     if (dd->ndim > 1)
2679     {
2680         dim0 = dd->dim[0];
2681         for (i = 0; i < dd->nc[dim0]; i++)
2682         {
2683             loc[dim0] = i;
2684             make_load_communicator(dd, 1, loc);
2685         }
2686     }
2687     if (dd->ndim > 2)
2688     {
2689         dim0 = dd->dim[0];
2690         for (i = 0; i < dd->nc[dim0]; i++)
2691         {
2692             loc[dim0] = i;
2693             dim1      = dd->dim[1];
2694             for (j = 0; j < dd->nc[dim1]; j++)
2695             {
2696                 loc[dim1] = j;
2697                 make_load_communicator(dd, 2, loc);
2698             }
2699         }
2700     }
2701
2702     if (debug)
2703     {
2704         fprintf(debug, "Finished making load communicators\n");
2705     }
2706 #endif
2707 }
2708
2709 /*! \brief Sets up the relation between neighboring domains and zones */
2710 static void setup_neighbor_relations(gmx_domdec_t *dd)
2711 {
2712     int                     d, dim, i, j, m;
2713     ivec                    tmp, s;
2714     gmx_domdec_zones_t     *zones;
2715     gmx_domdec_ns_ranges_t *izone;
2716
2717     for (d = 0; d < dd->ndim; d++)
2718     {
2719         dim = dd->dim[d];
2720         copy_ivec(dd->ci, tmp);
2721         tmp[dim]           = (tmp[dim] + 1) % dd->nc[dim];
2722         dd->neighbor[d][0] = ddcoord2ddnodeid(dd, tmp);
2723         copy_ivec(dd->ci, tmp);
2724         tmp[dim]           = (tmp[dim] - 1 + dd->nc[dim]) % dd->nc[dim];
2725         dd->neighbor[d][1] = ddcoord2ddnodeid(dd, tmp);
2726         if (debug)
2727         {
2728             fprintf(debug, "DD rank %d neighbor ranks in dir %d are + %d - %d\n",
2729                     dd->rank, dim,
2730                     dd->neighbor[d][0],
2731                     dd->neighbor[d][1]);
2732         }
2733     }
2734
2735     int nzone  = (1 << dd->ndim);
2736     int nizone = (1 << std::max(dd->ndim - 1, 0));
2737     assert(nizone >= 1 && nizone <= DD_MAXIZONE);
2738
2739     zones = &dd->comm->zones;
2740
2741     for (i = 0; i < nzone; i++)
2742     {
2743         m = 0;
2744         clear_ivec(zones->shift[i]);
2745         for (d = 0; d < dd->ndim; d++)
2746         {
2747             zones->shift[i][dd->dim[d]] = dd_zo[i][m++];
2748         }
2749     }
2750
2751     zones->n = nzone;
2752     for (i = 0; i < nzone; i++)
2753     {
2754         for (d = 0; d < DIM; d++)
2755         {
2756             s[d] = dd->ci[d] - zones->shift[i][d];
2757             if (s[d] < 0)
2758             {
2759                 s[d] += dd->nc[d];
2760             }
2761             else if (s[d] >= dd->nc[d])
2762             {
2763                 s[d] -= dd->nc[d];
2764             }
2765         }
2766     }
2767     zones->nizone = nizone;
2768     for (i = 0; i < zones->nizone; i++)
2769     {
2770         assert(ddNonbondedZonePairRanges[i][0] == i);
2771
2772         izone     = &zones->izone[i];
2773         /* dd_zp3 is for 3D decomposition, for fewer dimensions use only
2774          * j-zones up to nzone.
2775          */
2776         izone->j0 = std::min(ddNonbondedZonePairRanges[i][1], nzone);
2777         izone->j1 = std::min(ddNonbondedZonePairRanges[i][2], nzone);
2778         for (dim = 0; dim < DIM; dim++)
2779         {
2780             if (dd->nc[dim] == 1)
2781             {
2782                 /* All shifts should be allowed */
2783                 izone->shift0[dim] = -1;
2784                 izone->shift1[dim] = 1;
2785             }
2786             else
2787             {
2788                 /* Determine the min/max j-zone shift wrt the i-zone */
2789                 izone->shift0[dim] = 1;
2790                 izone->shift1[dim] = -1;
2791                 for (j = izone->j0; j < izone->j1; j++)
2792                 {
2793                     int shift_diff = zones->shift[j][dim] - zones->shift[i][dim];
2794                     if (shift_diff < izone->shift0[dim])
2795                     {
2796                         izone->shift0[dim] = shift_diff;
2797                     }
2798                     if (shift_diff > izone->shift1[dim])
2799                     {
2800                         izone->shift1[dim] = shift_diff;
2801                     }
2802                 }
2803             }
2804         }
2805     }
2806
2807     if (!isDlbDisabled(dd->comm))
2808     {
2809         dd->comm->cellsizesWithDlb.resize(dd->ndim);
2810     }
2811
2812     if (dd->comm->bRecordLoad)
2813     {
2814         make_load_communicators(dd);
2815     }
2816 }
2817
2818 static void make_pp_communicator(FILE                 *fplog,
2819                                  gmx_domdec_t         *dd,
2820                                  t_commrec gmx_unused *cr,
2821                                  int gmx_unused        reorder)
2822 {
2823 #if GMX_MPI
2824     gmx_domdec_comm_t *comm;
2825     int                rank, *buf;
2826     ivec               periods;
2827     MPI_Comm           comm_cart;
2828
2829     comm = dd->comm;
2830
2831     if (comm->bCartesianPP)
2832     {
2833         /* Set up cartesian communication for the particle-particle part */
2834         if (fplog)
2835         {
2836             fprintf(fplog, "Will use a Cartesian communicator: %d x %d x %d\n",
2837                     dd->nc[XX], dd->nc[YY], dd->nc[ZZ]);
2838         }
2839
2840         for (int i = 0; i < DIM; i++)
2841         {
2842             periods[i] = TRUE;
2843         }
2844         MPI_Cart_create(cr->mpi_comm_mygroup, DIM, dd->nc, periods, reorder,
2845                         &comm_cart);
2846         /* We overwrite the old communicator with the new cartesian one */
2847         cr->mpi_comm_mygroup = comm_cart;
2848     }
2849
2850     dd->mpi_comm_all = cr->mpi_comm_mygroup;
2851     MPI_Comm_rank(dd->mpi_comm_all, &dd->rank);
2852
2853     if (comm->bCartesianPP_PME)
2854     {
2855         /* Since we want to use the original cartesian setup for sim,
2856          * and not the one after split, we need to make an index.
2857          */
2858         snew(comm->ddindex2ddnodeid, dd->nnodes);
2859         comm->ddindex2ddnodeid[dd_index(dd->nc, dd->ci)] = dd->rank;
2860         gmx_sumi(dd->nnodes, comm->ddindex2ddnodeid, cr);
2861         /* Get the rank of the DD master,
2862          * above we made sure that the master node is a PP node.
2863          */
2864         if (MASTER(cr))
2865         {
2866             rank = dd->rank;
2867         }
2868         else
2869         {
2870             rank = 0;
2871         }
2872         MPI_Allreduce(&rank, &dd->masterrank, 1, MPI_INT, MPI_SUM, dd->mpi_comm_all);
2873     }
2874     else if (comm->bCartesianPP)
2875     {
2876         if (cr->npmenodes == 0)
2877         {
2878             /* The PP communicator is also
2879              * the communicator for this simulation
2880              */
2881             cr->mpi_comm_mysim = cr->mpi_comm_mygroup;
2882         }
2883         cr->nodeid = dd->rank;
2884
2885         MPI_Cart_coords(dd->mpi_comm_all, dd->rank, DIM, dd->ci);
2886
2887         /* We need to make an index to go from the coordinates
2888          * to the nodeid of this simulation.
2889          */
2890         snew(comm->ddindex2simnodeid, dd->nnodes);
2891         snew(buf, dd->nnodes);
2892         if (thisRankHasDuty(cr, DUTY_PP))
2893         {
2894             buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
2895         }
2896         /* Communicate the ddindex to simulation nodeid index */
2897         MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
2898                       cr->mpi_comm_mysim);
2899         sfree(buf);
2900
2901         /* Determine the master coordinates and rank.
2902          * The DD master should be the same node as the master of this sim.
2903          */
2904         for (int i = 0; i < dd->nnodes; i++)
2905         {
2906             if (comm->ddindex2simnodeid[i] == 0)
2907             {
2908                 ddindex2xyz(dd->nc, i, dd->master_ci);
2909                 MPI_Cart_rank(dd->mpi_comm_all, dd->master_ci, &dd->masterrank);
2910             }
2911         }
2912         if (debug)
2913         {
2914             fprintf(debug, "The master rank is %d\n", dd->masterrank);
2915         }
2916     }
2917     else
2918     {
2919         /* No Cartesian communicators */
2920         /* We use the rank in dd->comm->all as DD index */
2921         ddindex2xyz(dd->nc, dd->rank, dd->ci);
2922         /* The simulation master nodeid is 0, so the DD master rank is also 0 */
2923         dd->masterrank = 0;
2924         clear_ivec(dd->master_ci);
2925     }
2926 #endif
2927
2928     if (fplog)
2929     {
2930         fprintf(fplog,
2931                 "Domain decomposition rank %d, coordinates %d %d %d\n\n",
2932                 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
2933     }
2934     if (debug)
2935     {
2936         fprintf(debug,
2937                 "Domain decomposition rank %d, coordinates %d %d %d\n\n",
2938                 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
2939     }
2940 }
2941
2942 static void receive_ddindex2simnodeid(gmx_domdec_t         *dd,
2943                                       t_commrec            *cr)
2944 {
2945 #if GMX_MPI
2946     gmx_domdec_comm_t *comm = dd->comm;
2947
2948     if (!comm->bCartesianPP_PME && comm->bCartesianPP)
2949     {
2950         int *buf;
2951         snew(comm->ddindex2simnodeid, dd->nnodes);
2952         snew(buf, dd->nnodes);
2953         if (thisRankHasDuty(cr, DUTY_PP))
2954         {
2955             buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
2956         }
2957         /* Communicate the ddindex to simulation nodeid index */
2958         MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
2959                       cr->mpi_comm_mysim);
2960         sfree(buf);
2961     }
2962 #else
2963     GMX_UNUSED_VALUE(dd);
2964     GMX_UNUSED_VALUE(cr);
2965 #endif
2966 }
2967
2968 static void split_communicator(FILE *fplog, t_commrec *cr, gmx_domdec_t *dd,
2969                                DdRankOrder gmx_unused rankOrder,
2970                                int gmx_unused reorder)
2971 {
2972     gmx_domdec_comm_t *comm;
2973     int                i;
2974     gmx_bool           bDiv[DIM];
2975 #if GMX_MPI
2976     MPI_Comm           comm_cart;
2977 #endif
2978
2979     comm = dd->comm;
2980
2981     if (comm->bCartesianPP)
2982     {
2983         for (i = 1; i < DIM; i++)
2984         {
2985             bDiv[i] = ((cr->npmenodes*dd->nc[i]) % (dd->nnodes) == 0);
2986         }
2987         if (bDiv[YY] || bDiv[ZZ])
2988         {
2989             comm->bCartesianPP_PME = TRUE;
2990             /* If we have 2D PME decomposition, which is always in x+y,
2991              * we stack the PME only nodes in z.
2992              * Otherwise we choose the direction that provides the thinnest slab
2993              * of PME only nodes as this will have the least effect
2994              * on the PP communication.
2995              * But for the PME communication the opposite might be better.
2996              */
2997             if (bDiv[ZZ] && (comm->npmenodes_y > 1 ||
2998                              !bDiv[YY] ||
2999                              dd->nc[YY] > dd->nc[ZZ]))
3000             {
3001                 comm->cartpmedim = ZZ;
3002             }
3003             else
3004             {
3005                 comm->cartpmedim = YY;
3006             }
3007             comm->ntot[comm->cartpmedim]
3008                 += (cr->npmenodes*dd->nc[comm->cartpmedim])/dd->nnodes;
3009         }
3010         else if (fplog)
3011         {
3012             fprintf(fplog, "Number of PME-only ranks (%d) is not a multiple of nx*ny (%d*%d) or nx*nz (%d*%d)\n", cr->npmenodes, dd->nc[XX], dd->nc[YY], dd->nc[XX], dd->nc[ZZ]);
3013             fprintf(fplog,
3014                     "Will not use a Cartesian communicator for PP <-> PME\n\n");
3015         }
3016     }
3017
3018     if (comm->bCartesianPP_PME)
3019     {
3020 #if GMX_MPI
3021         int  rank;
3022         ivec periods;
3023
3024         if (fplog)
3025         {
3026             fprintf(fplog, "Will use a Cartesian communicator for PP <-> PME: %d x %d x %d\n", comm->ntot[XX], comm->ntot[YY], comm->ntot[ZZ]);
3027         }
3028
3029         for (i = 0; i < DIM; i++)
3030         {
3031             periods[i] = TRUE;
3032         }
3033         MPI_Cart_create(cr->mpi_comm_mysim, DIM, comm->ntot, periods, reorder,
3034                         &comm_cart);
3035         MPI_Comm_rank(comm_cart, &rank);
3036         if (MASTER(cr) && rank != 0)
3037         {
3038             gmx_fatal(FARGS, "MPI rank 0 was renumbered by MPI_Cart_create, we do not allow this");
3039         }
3040
3041         /* With this assigment we loose the link to the original communicator
3042          * which will usually be MPI_COMM_WORLD, unless have multisim.
3043          */
3044         cr->mpi_comm_mysim = comm_cart;
3045         cr->sim_nodeid     = rank;
3046
3047         MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, dd->ci);
3048
3049         if (fplog)
3050         {
3051             fprintf(fplog, "Cartesian rank %d, coordinates %d %d %d\n\n",
3052                     cr->sim_nodeid, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
3053         }
3054
3055         if (dd->ci[comm->cartpmedim] < dd->nc[comm->cartpmedim])
3056         {
3057             cr->duty = DUTY_PP;
3058         }
3059         if (cr->npmenodes == 0 ||
3060             dd->ci[comm->cartpmedim] >= dd->nc[comm->cartpmedim])
3061         {
3062             cr->duty = DUTY_PME;
3063         }
3064
3065         /* Split the sim communicator into PP and PME only nodes */
3066         MPI_Comm_split(cr->mpi_comm_mysim,
3067                        getThisRankDuties(cr),
3068                        dd_index(comm->ntot, dd->ci),
3069                        &cr->mpi_comm_mygroup);
3070 #endif
3071     }
3072     else
3073     {
3074         switch (rankOrder)
3075         {
3076             case DdRankOrder::pp_pme:
3077                 if (fplog)
3078                 {
3079                     fprintf(fplog, "Order of the ranks: PP first, PME last\n");
3080                 }
3081                 break;
3082             case DdRankOrder::interleave:
3083                 /* Interleave the PP-only and PME-only ranks */
3084                 if (fplog)
3085                 {
3086                     fprintf(fplog, "Interleaving PP and PME ranks\n");
3087                 }
3088                 comm->pmenodes = dd_interleaved_pme_ranks(dd);
3089                 break;
3090             case DdRankOrder::cartesian:
3091                 break;
3092             default:
3093                 gmx_fatal(FARGS, "Invalid ddRankOrder=%d", static_cast<int>(rankOrder));
3094         }
3095
3096         if (dd_simnode2pmenode(dd, cr, cr->sim_nodeid) == -1)
3097         {
3098             cr->duty = DUTY_PME;
3099         }
3100         else
3101         {
3102             cr->duty = DUTY_PP;
3103         }
3104 #if GMX_MPI
3105         /* Split the sim communicator into PP and PME only nodes */
3106         MPI_Comm_split(cr->mpi_comm_mysim,
3107                        getThisRankDuties(cr),
3108                        cr->nodeid,
3109                        &cr->mpi_comm_mygroup);
3110         MPI_Comm_rank(cr->mpi_comm_mygroup, &cr->nodeid);
3111 #endif
3112     }
3113
3114     if (fplog)
3115     {
3116         fprintf(fplog, "This rank does only %s work.\n\n",
3117                 thisRankHasDuty(cr, DUTY_PP) ? "particle-particle" : "PME-mesh");
3118     }
3119 }
3120
3121 /*! \brief Generates the MPI communicators for domain decomposition */
3122 static void make_dd_communicators(FILE *fplog, t_commrec *cr,
3123                                   gmx_domdec_t *dd, DdRankOrder ddRankOrder)
3124 {
3125     gmx_domdec_comm_t *comm;
3126     int                CartReorder;
3127
3128     comm = dd->comm;
3129
3130     copy_ivec(dd->nc, comm->ntot);
3131
3132     comm->bCartesianPP     = (ddRankOrder == DdRankOrder::cartesian);
3133     comm->bCartesianPP_PME = FALSE;
3134
3135     /* Reorder the nodes by default. This might change the MPI ranks.
3136      * Real reordering is only supported on very few architectures,
3137      * Blue Gene is one of them.
3138      */
3139     CartReorder = (getenv("GMX_NO_CART_REORDER") == nullptr);
3140
3141     if (cr->npmenodes > 0)
3142     {
3143         /* Split the communicator into a PP and PME part */
3144         split_communicator(fplog, cr, dd, ddRankOrder, CartReorder);
3145         if (comm->bCartesianPP_PME)
3146         {
3147             /* We (possibly) reordered the nodes in split_communicator,
3148              * so it is no longer required in make_pp_communicator.
3149              */
3150             CartReorder = FALSE;
3151         }
3152     }
3153     else
3154     {
3155         /* All nodes do PP and PME */
3156 #if GMX_MPI
3157         /* We do not require separate communicators */
3158         cr->mpi_comm_mygroup = cr->mpi_comm_mysim;
3159 #endif
3160     }
3161
3162     if (thisRankHasDuty(cr, DUTY_PP))
3163     {
3164         /* Copy or make a new PP communicator */
3165         make_pp_communicator(fplog, dd, cr, CartReorder);
3166     }
3167     else
3168     {
3169         receive_ddindex2simnodeid(dd, cr);
3170     }
3171
3172     if (!thisRankHasDuty(cr, DUTY_PME))
3173     {
3174         /* Set up the commnuication to our PME node */
3175         dd->pme_nodeid           = dd_simnode2pmenode(dd, cr, cr->sim_nodeid);
3176         dd->pme_receive_vir_ener = receive_vir_ener(dd, cr);
3177         if (debug)
3178         {
3179             fprintf(debug, "My pme_nodeid %d receive ener %s\n",
3180                     dd->pme_nodeid, gmx::boolToString(dd->pme_receive_vir_ener));
3181         }
3182     }
3183     else
3184     {
3185         dd->pme_nodeid = -1;
3186     }
3187
3188     if (DDMASTER(dd))
3189     {
3190         dd->ma = gmx::compat::make_unique<AtomDistribution>(dd->nc,
3191                                                             comm->cgs_gl.nr,
3192                                                             comm->cgs_gl.index[comm->cgs_gl.nr]);
3193     }
3194 }
3195
3196 static real *get_slb_frac(FILE *fplog, const char *dir, int nc, const char *size_string)
3197 {
3198     real  *slb_frac, tot;
3199     int    i, n;
3200     double dbl;
3201
3202     slb_frac = nullptr;
3203     if (nc > 1 && size_string != nullptr)
3204     {
3205         if (fplog)
3206         {
3207             fprintf(fplog, "Using static load balancing for the %s direction\n",
3208                     dir);
3209         }
3210         snew(slb_frac, nc);
3211         tot = 0;
3212         for (i = 0; i < nc; i++)
3213         {
3214             dbl = 0;
3215             sscanf(size_string, "%20lf%n", &dbl, &n);
3216             if (dbl == 0)
3217             {
3218                 gmx_fatal(FARGS, "Incorrect or not enough DD cell size entries for direction %s: '%s'", dir, size_string);
3219             }
3220             slb_frac[i]  = dbl;
3221             size_string += n;
3222             tot         += slb_frac[i];
3223         }
3224         /* Normalize */
3225         if (fplog)
3226         {
3227             fprintf(fplog, "Relative cell sizes:");
3228         }
3229         for (i = 0; i < nc; i++)
3230         {
3231             slb_frac[i] /= tot;
3232             if (fplog)
3233             {
3234                 fprintf(fplog, " %5.3f", slb_frac[i]);
3235             }
3236         }
3237         if (fplog)
3238         {
3239             fprintf(fplog, "\n");
3240         }
3241     }
3242
3243     return slb_frac;
3244 }
3245
3246 static int multi_body_bondeds_count(const gmx_mtop_t *mtop)
3247 {
3248     int                  n, nmol, ftype;
3249     gmx_mtop_ilistloop_t iloop;
3250     const t_ilist       *il;
3251
3252     n     = 0;
3253     iloop = gmx_mtop_ilistloop_init(mtop);
3254     while (gmx_mtop_ilistloop_next(iloop, &il, &nmol))
3255     {
3256         for (ftype = 0; ftype < F_NRE; ftype++)
3257         {
3258             if ((interaction_function[ftype].flags & IF_BOND) &&
3259                 NRAL(ftype) >  2)
3260             {
3261                 n += nmol*il[ftype].nr/(1 + NRAL(ftype));
3262             }
3263         }
3264     }
3265
3266     return n;
3267 }
3268
3269 static int dd_getenv(FILE *fplog, const char *env_var, int def)
3270 {
3271     char *val;
3272     int   nst;
3273
3274     nst = def;
3275     val = getenv(env_var);
3276     if (val)
3277     {
3278         if (sscanf(val, "%20d", &nst) <= 0)
3279         {
3280             nst = 1;
3281         }
3282         if (fplog)
3283         {
3284             fprintf(fplog, "Found env.var. %s = %s, using value %d\n",
3285                     env_var, val, nst);
3286         }
3287     }
3288
3289     return nst;
3290 }
3291
3292 static void dd_warning(const t_commrec *cr, FILE *fplog, const char *warn_string)
3293 {
3294     if (MASTER(cr))
3295     {
3296         fprintf(stderr, "\n%s\n", warn_string);
3297     }
3298     if (fplog)
3299     {
3300         fprintf(fplog, "\n%s\n", warn_string);
3301     }
3302 }
3303
3304 static void check_dd_restrictions(t_commrec *cr, const gmx_domdec_t *dd,
3305                                   const t_inputrec *ir, FILE *fplog)
3306 {
3307     if (ir->ePBC == epbcSCREW &&
3308         (dd->nc[XX] == 1 || dd->nc[YY] > 1 || dd->nc[ZZ] > 1))
3309     {
3310         gmx_fatal(FARGS, "With pbc=%s can only do domain decomposition in the x-direction", epbc_names[ir->ePBC]);
3311     }
3312
3313     if (ir->ns_type == ensSIMPLE)
3314     {
3315         gmx_fatal(FARGS, "Domain decomposition does not support simple neighbor searching, use grid searching or run with one MPI rank");
3316     }
3317
3318     if (ir->nstlist == 0)
3319     {
3320         gmx_fatal(FARGS, "Domain decomposition does not work with nstlist=0");
3321     }
3322
3323     if (ir->comm_mode == ecmANGULAR && ir->ePBC != epbcNONE)
3324     {
3325         dd_warning(cr, fplog, "comm-mode angular will give incorrect results when the comm group partially crosses a periodic boundary");
3326     }
3327 }
3328
3329 static real average_cellsize_min(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
3330 {
3331     int  di, d;
3332     real r;
3333
3334     r = ddbox->box_size[XX];
3335     for (di = 0; di < dd->ndim; di++)
3336     {
3337         d = dd->dim[di];
3338         /* Check using the initial average cell size */
3339         r = std::min(r, ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
3340     }
3341
3342     return r;
3343 }
3344
3345 /*! \brief Depending on the DLB initial value return the DLB switched off state or issue an error.
3346  */
3347 static int forceDlbOffOrBail(int                cmdlineDlbState,
3348                              const std::string &reasonStr,
3349                              t_commrec         *cr,
3350                              FILE              *fplog)
3351 {
3352     std::string dlbNotSupportedErr  = "Dynamic load balancing requested, but ";
3353     std::string dlbDisableNote      = "NOTE: disabling dynamic load balancing as ";
3354
3355     if (cmdlineDlbState == edlbsOnUser)
3356     {
3357         gmx_fatal(FARGS, "%s", (dlbNotSupportedErr + reasonStr).c_str());
3358     }
3359     else if (cmdlineDlbState == edlbsOffCanTurnOn)
3360     {
3361         dd_warning(cr, fplog, (dlbDisableNote + reasonStr + "\n").c_str());
3362     }
3363     return edlbsOffForever;
3364 }
3365
3366 /*! \brief Return the dynamic load balancer's initial state based on initial conditions and user inputs.
3367  *
3368  * This function parses the parameters of "-dlb" command line option setting
3369  * corresponding state values. Then it checks the consistency of the determined
3370  * state with other run parameters and settings. As a result, the initial state
3371  * may be altered or an error may be thrown if incompatibility of options is detected.
3372  *
3373  * \param [in] fplog       Pointer to mdrun log file.
3374  * \param [in] cr          Pointer to MPI communication object.
3375  * \param [in] dlbOption   Enum value for the DLB option.
3376  * \param [in] bRecordLoad True if the load balancer is recording load information.
3377  * \param [in] mdrunOptions  Options for mdrun.
3378  * \param [in] ir          Pointer mdrun to input parameters.
3379  * \returns                DLB initial/startup state.
3380  */
3381 static int determineInitialDlbState(FILE *fplog, t_commrec *cr,
3382                                     DlbOption dlbOption, gmx_bool bRecordLoad,
3383                                     const MdrunOptions &mdrunOptions,
3384                                     const t_inputrec *ir)
3385 {
3386     int dlbState = edlbsOffCanTurnOn;
3387
3388     switch (dlbOption)
3389     {
3390         case DlbOption::turnOnWhenUseful: dlbState = edlbsOffCanTurnOn; break;
3391         case DlbOption::no:               dlbState = edlbsOffUser;      break;
3392         case DlbOption::yes:              dlbState = edlbsOnUser;       break;
3393         default: gmx_incons("Invalid dlbOption enum value");
3394     }
3395
3396     /* Reruns don't support DLB: bail or override auto mode */
3397     if (mdrunOptions.rerun)
3398     {
3399         std::string reasonStr = "it is not supported in reruns.";
3400         return forceDlbOffOrBail(dlbState, reasonStr, cr, fplog);
3401     }
3402
3403     /* Unsupported integrators */
3404     if (!EI_DYNAMICS(ir->eI))
3405     {
3406         auto reasonStr = gmx::formatString("it is only supported with dynamics, not with integrator '%s'.", EI(ir->eI));
3407         return forceDlbOffOrBail(dlbState, reasonStr, cr, fplog);
3408     }
3409
3410     /* Without cycle counters we can't time work to balance on */
3411     if (!bRecordLoad)
3412     {
3413         std::string reasonStr = "cycle counters unsupported or not enabled in the operating system kernel.";
3414         return forceDlbOffOrBail(dlbState, reasonStr, cr, fplog);
3415     }
3416
3417     if (mdrunOptions.reproducible)
3418     {
3419         std::string reasonStr = "you started a reproducible run.";
3420         switch (dlbState)
3421         {
3422             case edlbsOffUser:
3423                 break;
3424             case edlbsOffForever:
3425                 GMX_RELEASE_ASSERT(false, "edlbsOffForever is not a valid initial state");
3426                 break;
3427             case edlbsOffCanTurnOn:
3428                 return forceDlbOffOrBail(dlbState, reasonStr, cr, fplog);
3429             case edlbsOnCanTurnOff:
3430                 GMX_RELEASE_ASSERT(false, "edlbsOffCanTurnOff is not a valid initial state");
3431                 break;
3432             case edlbsOnUser:
3433                 return forceDlbOffOrBail(dlbState, reasonStr + " In load balanced runs binary reproducibility cannot be ensured.", cr, fplog);
3434             default:
3435                 gmx_fatal(FARGS, "Death horror: undefined case (%d) for load balancing choice", dlbState);
3436         }
3437     }
3438
3439     return dlbState;
3440 }
3441
3442 static void set_dd_dim(FILE *fplog, gmx_domdec_t *dd)
3443 {
3444     dd->ndim = 0;
3445     if (getenv("GMX_DD_ORDER_ZYX") != nullptr)
3446     {
3447         /* Decomposition order z,y,x */
3448         if (fplog)
3449         {
3450             fprintf(fplog, "Using domain decomposition order z, y, x\n");
3451         }
3452         for (int dim = DIM-1; dim >= 0; dim--)
3453         {
3454             if (dd->nc[dim] > 1)
3455             {
3456                 dd->dim[dd->ndim++] = dim;
3457             }
3458         }
3459     }
3460     else
3461     {
3462         /* Decomposition order x,y,z */
3463         for (int dim = 0; dim < DIM; dim++)
3464         {
3465             if (dd->nc[dim] > 1)
3466             {
3467                 dd->dim[dd->ndim++] = dim;
3468             }
3469         }
3470     }
3471
3472     if (dd->ndim == 0)
3473     {
3474         /* Set dim[0] to avoid extra checks on ndim in several places */
3475         dd->dim[0] = XX;
3476     }
3477 }
3478
3479 static gmx_domdec_comm_t *init_dd_comm()
3480 {
3481     gmx_domdec_comm_t *comm = new gmx_domdec_comm_t;
3482
3483     comm->n_load_have      = 0;
3484     comm->n_load_collect   = 0;
3485
3486     comm->haveTurnedOffDlb = false;
3487
3488     for (int i = 0; i < static_cast<int>(DDAtomRanges::Type::Number); i++)
3489     {
3490         comm->sum_nat[i] = 0;
3491     }
3492     comm->ndecomp   = 0;
3493     comm->nload     = 0;
3494     comm->load_step = 0;
3495     comm->load_sum  = 0;
3496     comm->load_max  = 0;
3497     clear_ivec(comm->load_lim);
3498     comm->load_mdf  = 0;
3499     comm->load_pme  = 0;
3500
3501     /* This should be replaced by a unique pointer */
3502     comm->balanceRegion = ddBalanceRegionAllocate();
3503
3504     return comm;
3505 }
3506
3507 /*! \brief Set the cell size and interaction limits, as well as the DD grid */
3508 static void set_dd_limits_and_grid(FILE *fplog, t_commrec *cr, gmx_domdec_t *dd,
3509                                    const DomdecOptions &options,
3510                                    const MdrunOptions &mdrunOptions,
3511                                    const gmx_mtop_t *mtop,
3512                                    const t_inputrec *ir,
3513                                    const matrix box,
3514                                    gmx::ArrayRef<const gmx::RVec> xGlobal,
3515                                    gmx_ddbox_t *ddbox)
3516 {
3517     real               r_bonded         = -1;
3518     real               r_bonded_limit   = -1;
3519     const real         tenPercentMargin = 1.1;
3520     gmx_domdec_comm_t *comm             = dd->comm;
3521
3522     dd->npbcdim   = ePBC2npbcdim(ir->ePBC);
3523     dd->bScrewPBC = (ir->ePBC == epbcSCREW);
3524
3525     dd->pme_recv_f_alloc = 0;
3526     dd->pme_recv_f_buf   = nullptr;
3527
3528     /* Initialize to GPU share count to 0, might change later */
3529     comm->nrank_gpu_shared = 0;
3530
3531     comm->dlbState         = determineInitialDlbState(fplog, cr, options.dlbOption, comm->bRecordLoad, mdrunOptions, ir);
3532     dd_dlb_set_should_check_whether_to_turn_dlb_on(dd, TRUE);
3533     /* To consider turning DLB on after 2*nstlist steps we need to check
3534      * at partitioning count 3. Thus we need to increase the first count by 2.
3535      */
3536     comm->ddPartioningCountFirstDlbOff += 2;
3537
3538     if (fplog)
3539     {
3540         fprintf(fplog, "Dynamic load balancing: %s\n",
3541                 edlbs_names[comm->dlbState]);
3542     }
3543     comm->bPMELoadBalDLBLimits = FALSE;
3544
3545     /* Allocate the charge group/atom sorting struct */
3546     comm->sort = gmx::compat::make_unique<gmx_domdec_sort_t>();
3547
3548     comm->bCGs = (ncg_mtop(mtop) < mtop->natoms);
3549
3550     comm->bInterCGBondeds = ((ncg_mtop(mtop) > gmx_mtop_num_molecules(*mtop)) ||
3551                              mtop->bIntermolecularInteractions);
3552     if (comm->bInterCGBondeds)
3553     {
3554         comm->bInterCGMultiBody = (multi_body_bondeds_count(mtop) > 0);
3555     }
3556     else
3557     {
3558         comm->bInterCGMultiBody = FALSE;
3559     }
3560
3561     dd->bInterCGcons    = gmx::inter_charge_group_constraints(*mtop);
3562     dd->bInterCGsettles = gmx::inter_charge_group_settles(*mtop);
3563
3564     if (ir->rlist == 0)
3565     {
3566         /* Set the cut-off to some very large value,
3567          * so we don't need if statements everywhere in the code.
3568          * We use sqrt, since the cut-off is squared in some places.
3569          */
3570         comm->cutoff   = GMX_CUTOFF_INF;
3571     }
3572     else
3573     {
3574         comm->cutoff   = ir->rlist;
3575     }
3576     comm->cutoff_mbody = 0;
3577
3578     comm->cellsize_limit = 0;
3579     comm->bBondComm      = FALSE;
3580
3581     /* Atoms should be able to move by up to half the list buffer size (if > 0)
3582      * within nstlist steps. Since boundaries are allowed to displace by half
3583      * a cell size, DD cells should be at least the size of the list buffer.
3584      */
3585     comm->cellsize_limit = std::max(comm->cellsize_limit,
3586                                     ir->rlist - std::max(ir->rvdw, ir->rcoulomb));
3587
3588     if (comm->bInterCGBondeds)
3589     {
3590         if (options.minimumCommunicationRange > 0)
3591         {
3592             comm->cutoff_mbody = options.minimumCommunicationRange;
3593             if (options.useBondedCommunication)
3594             {
3595                 comm->bBondComm = (comm->cutoff_mbody > comm->cutoff);
3596             }
3597             else
3598             {
3599                 comm->cutoff = std::max(comm->cutoff, comm->cutoff_mbody);
3600             }
3601             r_bonded_limit = comm->cutoff_mbody;
3602         }
3603         else if (ir->bPeriodicMols)
3604         {
3605             /* Can not easily determine the required cut-off */
3606             dd_warning(cr, fplog, "NOTE: Periodic molecules are present in this system. Because of this, the domain decomposition algorithm cannot easily determine the minimum cell size that it requires for treating bonded interactions. Instead, domain decomposition will assume that half the non-bonded cut-off will be a suitable lower bound.\n");
3607             comm->cutoff_mbody = comm->cutoff/2;
3608             r_bonded_limit     = comm->cutoff_mbody;
3609         }
3610         else
3611         {
3612             real r_2b, r_mb;
3613
3614             if (MASTER(cr))
3615             {
3616                 dd_bonded_cg_distance(fplog, mtop, ir, as_rvec_array(xGlobal.data()), box,
3617                                       options.checkBondedInteractions,
3618                                       &r_2b, &r_mb);
3619             }
3620             gmx_bcast(sizeof(r_2b), &r_2b, cr);
3621             gmx_bcast(sizeof(r_mb), &r_mb, cr);
3622
3623             /* We use an initial margin of 10% for the minimum cell size,
3624              * except when we are just below the non-bonded cut-off.
3625              */
3626             if (options.useBondedCommunication)
3627             {
3628                 if (std::max(r_2b, r_mb) > comm->cutoff)
3629                 {
3630                     r_bonded        = std::max(r_2b, r_mb);
3631                     r_bonded_limit  = tenPercentMargin*r_bonded;
3632                     comm->bBondComm = TRUE;
3633                 }
3634                 else
3635                 {
3636                     r_bonded       = r_mb;
3637                     r_bonded_limit = std::min(tenPercentMargin*r_bonded, comm->cutoff);
3638                 }
3639                 /* We determine cutoff_mbody later */
3640             }
3641             else
3642             {
3643                 /* No special bonded communication,
3644                  * simply increase the DD cut-off.
3645                  */
3646                 r_bonded_limit     = tenPercentMargin*std::max(r_2b, r_mb);
3647                 comm->cutoff_mbody = r_bonded_limit;
3648                 comm->cutoff       = std::max(comm->cutoff, comm->cutoff_mbody);
3649             }
3650         }
3651         if (fplog)
3652         {
3653             fprintf(fplog,
3654                     "Minimum cell size due to bonded interactions: %.3f nm\n",
3655                     r_bonded_limit);
3656         }
3657         comm->cellsize_limit = std::max(comm->cellsize_limit, r_bonded_limit);
3658     }
3659
3660     real rconstr = 0;
3661     if (dd->bInterCGcons && options.constraintCommunicationRange <= 0)
3662     {
3663         /* There is a cell size limit due to the constraints (P-LINCS) */
3664         rconstr = gmx::constr_r_max(fplog, mtop, ir);
3665         if (fplog)
3666         {
3667             fprintf(fplog,
3668                     "Estimated maximum distance required for P-LINCS: %.3f nm\n",
3669                     rconstr);
3670             if (rconstr > comm->cellsize_limit)
3671             {
3672                 fprintf(fplog, "This distance will limit the DD cell size, you can override this with -rcon\n");
3673             }
3674         }
3675     }
3676     else if (options.constraintCommunicationRange > 0 && fplog)
3677     {
3678         /* Here we do not check for dd->bInterCGcons,
3679          * because one can also set a cell size limit for virtual sites only
3680          * and at this point we don't know yet if there are intercg v-sites.
3681          */
3682         fprintf(fplog,
3683                 "User supplied maximum distance required for P-LINCS: %.3f nm\n",
3684                 options.constraintCommunicationRange);
3685         rconstr = options.constraintCommunicationRange;
3686     }
3687     comm->cellsize_limit = std::max(comm->cellsize_limit, rconstr);
3688
3689     comm->cgs_gl = gmx_mtop_global_cgs(mtop);
3690
3691     if (options.numCells[XX] > 0)
3692     {
3693         copy_ivec(options.numCells, dd->nc);
3694         set_dd_dim(fplog, dd);
3695         set_ddbox_cr(cr, &dd->nc, ir, box, xGlobal, ddbox);
3696
3697         if (options.numPmeRanks >= 0)
3698         {
3699             cr->npmenodes = options.numPmeRanks;
3700         }
3701         else
3702         {
3703             /* When the DD grid is set explicitly and -npme is set to auto,
3704              * don't use PME ranks. We check later if the DD grid is
3705              * compatible with the total number of ranks.
3706              */
3707             cr->npmenodes = 0;
3708         }
3709
3710         real acs = average_cellsize_min(dd, ddbox);
3711         if (acs < comm->cellsize_limit)
3712         {
3713             if (fplog)
3714             {
3715                 fprintf(fplog, "ERROR: The initial cell size (%f) is smaller than the cell size limit (%f)\n", acs, comm->cellsize_limit);
3716             }
3717             gmx_fatal_collective(FARGS, cr->mpi_comm_mysim, MASTER(cr),
3718                                  "The initial cell size (%f) is smaller than the cell size limit (%f), change options -dd, -rdd or -rcon, see the log file for details",
3719                                  acs, comm->cellsize_limit);
3720         }
3721     }
3722     else
3723     {
3724         set_ddbox_cr(cr, nullptr, ir, box, xGlobal, ddbox);
3725
3726         /* We need to choose the optimal DD grid and possibly PME nodes */
3727         real limit =
3728             dd_choose_grid(fplog, cr, dd, ir, mtop, box, ddbox,
3729                            options.numPmeRanks,
3730                            !isDlbDisabled(comm),
3731                            options.dlbScaling,
3732                            comm->cellsize_limit, comm->cutoff,
3733                            comm->bInterCGBondeds);
3734
3735         if (dd->nc[XX] == 0)
3736         {
3737             char     buf[STRLEN];
3738             gmx_bool bC = (dd->bInterCGcons && rconstr > r_bonded_limit);
3739             sprintf(buf, "Change the number of ranks or mdrun option %s%s%s",
3740                     !bC ? "-rdd" : "-rcon",
3741                     comm->dlbState != edlbsOffUser ? " or -dds" : "",
3742                     bC ? " or your LINCS settings" : "");
3743
3744             gmx_fatal_collective(FARGS, cr->mpi_comm_mysim, MASTER(cr),
3745                                  "There is no domain decomposition for %d ranks that is compatible with the given box and a minimum cell size of %g nm\n"
3746                                  "%s\n"
3747                                  "Look in the log file for details on the domain decomposition",
3748                                  cr->nnodes-cr->npmenodes, limit, buf);
3749         }
3750         set_dd_dim(fplog, dd);
3751     }
3752
3753     if (fplog)
3754     {
3755         fprintf(fplog,
3756                 "Domain decomposition grid %d x %d x %d, separate PME ranks %d\n",
3757                 dd->nc[XX], dd->nc[YY], dd->nc[ZZ], cr->npmenodes);
3758     }
3759
3760     dd->nnodes = dd->nc[XX]*dd->nc[YY]*dd->nc[ZZ];
3761     if (cr->nnodes - dd->nnodes != cr->npmenodes)
3762     {
3763         gmx_fatal_collective(FARGS, cr->mpi_comm_mysim, MASTER(cr),
3764                              "The size of the domain decomposition grid (%d) does not match the number of ranks (%d). The total number of ranks is %d",
3765                              dd->nnodes, cr->nnodes - cr->npmenodes, cr->nnodes);
3766     }
3767     if (cr->npmenodes > dd->nnodes)
3768     {
3769         gmx_fatal_collective(FARGS, cr->mpi_comm_mysim, MASTER(cr),
3770                              "The number of separate PME ranks (%d) is larger than the number of PP ranks (%d), this is not supported.", cr->npmenodes, dd->nnodes);
3771     }
3772     if (cr->npmenodes > 0)
3773     {
3774         comm->npmenodes = cr->npmenodes;
3775     }
3776     else
3777     {
3778         comm->npmenodes = dd->nnodes;
3779     }
3780
3781     if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
3782     {
3783         /* The following choices should match those
3784          * in comm_cost_est in domdec_setup.c.
3785          * Note that here the checks have to take into account
3786          * that the decomposition might occur in a different order than xyz
3787          * (for instance through the env.var. GMX_DD_ORDER_ZYX),
3788          * in which case they will not match those in comm_cost_est,
3789          * but since that is mainly for testing purposes that's fine.
3790          */
3791         if (dd->ndim >= 2 && dd->dim[0] == XX && dd->dim[1] == YY &&
3792             comm->npmenodes > dd->nc[XX] && comm->npmenodes % dd->nc[XX] == 0 &&
3793             getenv("GMX_PMEONEDD") == nullptr)
3794         {
3795             comm->npmedecompdim = 2;
3796             comm->npmenodes_x   = dd->nc[XX];
3797             comm->npmenodes_y   = comm->npmenodes/comm->npmenodes_x;
3798         }
3799         else
3800         {
3801             /* In case nc is 1 in both x and y we could still choose to
3802              * decompose pme in y instead of x, but we use x for simplicity.
3803              */
3804             comm->npmedecompdim = 1;
3805             if (dd->dim[0] == YY)
3806             {
3807                 comm->npmenodes_x = 1;
3808                 comm->npmenodes_y = comm->npmenodes;
3809             }
3810             else
3811             {
3812                 comm->npmenodes_x = comm->npmenodes;
3813                 comm->npmenodes_y = 1;
3814             }
3815         }
3816         if (fplog)
3817         {
3818             fprintf(fplog, "PME domain decomposition: %d x %d x %d\n",
3819                     comm->npmenodes_x, comm->npmenodes_y, 1);
3820         }
3821     }
3822     else
3823     {
3824         comm->npmedecompdim = 0;
3825         comm->npmenodes_x   = 0;
3826         comm->npmenodes_y   = 0;
3827     }
3828
3829     snew(comm->slb_frac, DIM);
3830     if (isDlbDisabled(comm))
3831     {
3832         comm->slb_frac[XX] = get_slb_frac(fplog, "x", dd->nc[XX], options.cellSizeX);
3833         comm->slb_frac[YY] = get_slb_frac(fplog, "y", dd->nc[YY], options.cellSizeY);
3834         comm->slb_frac[ZZ] = get_slb_frac(fplog, "z", dd->nc[ZZ], options.cellSizeZ);
3835     }
3836
3837     if (comm->bInterCGBondeds && comm->cutoff_mbody == 0)
3838     {
3839         if (comm->bBondComm || !isDlbDisabled(comm))
3840         {
3841             /* Set the bonded communication distance to halfway
3842              * the minimum and the maximum,
3843              * since the extra communication cost is nearly zero.
3844              */
3845             real acs           = average_cellsize_min(dd, ddbox);
3846             comm->cutoff_mbody = 0.5*(r_bonded + acs);
3847             if (!isDlbDisabled(comm))
3848             {
3849                 /* Check if this does not limit the scaling */
3850                 comm->cutoff_mbody = std::min(comm->cutoff_mbody,
3851                                               options.dlbScaling*acs);
3852             }
3853             if (!comm->bBondComm)
3854             {
3855                 /* Without bBondComm do not go beyond the n.b. cut-off */
3856                 comm->cutoff_mbody = std::min(comm->cutoff_mbody, comm->cutoff);
3857                 if (comm->cellsize_limit >= comm->cutoff)
3858                 {
3859                     /* We don't loose a lot of efficieny
3860                      * when increasing it to the n.b. cut-off.
3861                      * It can even be slightly faster, because we need
3862                      * less checks for the communication setup.
3863                      */
3864                     comm->cutoff_mbody = comm->cutoff;
3865                 }
3866             }
3867             /* Check if we did not end up below our original limit */
3868             comm->cutoff_mbody = std::max(comm->cutoff_mbody, r_bonded_limit);
3869
3870             if (comm->cutoff_mbody > comm->cellsize_limit)
3871             {
3872                 comm->cellsize_limit = comm->cutoff_mbody;
3873             }
3874         }
3875         /* Without DLB and cutoff_mbody<cutoff, cutoff_mbody is dynamic */
3876     }
3877
3878     if (debug)
3879     {
3880         fprintf(debug, "Bonded atom communication beyond the cut-off: %s\n"
3881                 "cellsize limit %f\n",
3882                 gmx::boolToString(comm->bBondComm), comm->cellsize_limit);
3883     }
3884
3885     if (MASTER(cr))
3886     {
3887         check_dd_restrictions(cr, dd, ir, fplog);
3888     }
3889 }
3890
3891 static void set_dlb_limits(gmx_domdec_t *dd)
3892
3893 {
3894     for (int d = 0; d < dd->ndim; d++)
3895     {
3896         /* Set the number of pulses to the value for DLB */
3897         dd->comm->cd[d].ind.resize(dd->comm->cd[d].np_dlb);
3898
3899         dd->comm->cellsize_min[dd->dim[d]] =
3900             dd->comm->cellsize_min_dlb[dd->dim[d]];
3901     }
3902 }
3903
3904
3905 static void turn_on_dlb(FILE *fplog, const t_commrec *cr, int64_t step)
3906 {
3907     gmx_domdec_t      *dd           = cr->dd;
3908     gmx_domdec_comm_t *comm         = dd->comm;
3909
3910     real               cellsize_min = comm->cellsize_min[dd->dim[0]];
3911     for (int d = 1; d < dd->ndim; d++)
3912     {
3913         cellsize_min = std::min(cellsize_min, comm->cellsize_min[dd->dim[d]]);
3914     }
3915
3916     /* Turn off DLB if we're too close to the cell size limit. */
3917     if (cellsize_min < comm->cellsize_limit*1.05)
3918     {
3919         auto str = gmx::formatString("step %" PRId64 " Measured %.1f %% performance loss due to load imbalance, "
3920                                      "but the minimum cell size is smaller than 1.05 times the cell size limit."
3921                                      "Will no longer try dynamic load balancing.\n", step, dd_force_imb_perf_loss(dd)*100);
3922         dd_warning(cr, fplog, str.c_str());
3923
3924         comm->dlbState = edlbsOffForever;
3925         return;
3926     }
3927
3928     char buf[STRLEN];
3929     sprintf(buf, "step %" PRId64 " Turning on dynamic load balancing, because the performance loss due to load imbalance is %.1f %%.\n", step, dd_force_imb_perf_loss(dd)*100);
3930     dd_warning(cr, fplog, buf);
3931     comm->dlbState = edlbsOnCanTurnOff;
3932
3933     /* Store the non-DLB performance, so we can check if DLB actually
3934      * improves performance.
3935      */
3936     GMX_RELEASE_ASSERT(comm->cycl_n[ddCyclStep] > 0, "When we turned on DLB, we should have measured cycles");
3937     comm->cyclesPerStepBeforeDLB = comm->cycl[ddCyclStep]/comm->cycl_n[ddCyclStep];
3938
3939     set_dlb_limits(dd);
3940
3941     /* We can set the required cell size info here,
3942      * so we do not need to communicate this.
3943      * The grid is completely uniform.
3944      */
3945     for (int d = 0; d < dd->ndim; d++)
3946     {
3947         RowMaster *rowMaster = comm->cellsizesWithDlb[d].rowMaster.get();
3948
3949         if (rowMaster)
3950         {
3951             comm->load[d].sum_m = comm->load[d].sum;
3952
3953             int nc = dd->nc[dd->dim[d]];
3954             for (int i = 0; i < nc; i++)
3955             {
3956                 rowMaster->cellFrac[i] = i/static_cast<real>(nc);
3957                 if (d > 0)
3958                 {
3959                     rowMaster->bounds[i].cellFracLowerMax =  i     /static_cast<real>(nc);
3960                     rowMaster->bounds[i].cellFracUpperMin = (i + 1)/static_cast<real>(nc);
3961                 }
3962             }
3963             rowMaster->cellFrac[nc] = 1.0;
3964         }
3965     }
3966 }
3967
3968 static void turn_off_dlb(FILE *fplog, const t_commrec *cr, int64_t step)
3969 {
3970     gmx_domdec_t *dd = cr->dd;
3971
3972     char          buf[STRLEN];
3973     sprintf(buf, "step %" PRId64 " Turning off dynamic load balancing, because it is degrading performance.\n", step);
3974     dd_warning(cr, fplog, buf);
3975     dd->comm->dlbState                     = edlbsOffCanTurnOn;
3976     dd->comm->haveTurnedOffDlb             = true;
3977     dd->comm->ddPartioningCountFirstDlbOff = dd->ddp_count;
3978 }
3979
3980 static void turn_off_dlb_forever(FILE *fplog, const t_commrec *cr, int64_t step)
3981 {
3982     GMX_RELEASE_ASSERT(cr->dd->comm->dlbState == edlbsOffCanTurnOn, "Can only turn off DLB forever when it was in the can-turn-on state");
3983     char buf[STRLEN];
3984     sprintf(buf, "step %" PRId64 " Will no longer try dynamic load balancing, as it degraded performance.\n", step);
3985     dd_warning(cr, fplog, buf);
3986     cr->dd->comm->dlbState = edlbsOffForever;
3987 }
3988
3989 static char *init_bLocalCG(const gmx_mtop_t *mtop)
3990 {
3991     int   ncg, cg;
3992     char *bLocalCG;
3993
3994     ncg = ncg_mtop(mtop);
3995     snew(bLocalCG, ncg);
3996     for (cg = 0; cg < ncg; cg++)
3997     {
3998         bLocalCG[cg] = FALSE;
3999     }
4000
4001     return bLocalCG;
4002 }
4003
4004 void dd_init_bondeds(FILE *fplog,
4005                      gmx_domdec_t *dd,
4006                      const gmx_mtop_t *mtop,
4007                      const gmx_vsite_t *vsite,
4008                      const t_inputrec *ir,
4009                      gmx_bool bBCheck, cginfo_mb_t *cginfo_mb)
4010 {
4011     gmx_domdec_comm_t *comm;
4012
4013     dd_make_reverse_top(fplog, dd, mtop, vsite, ir, bBCheck);
4014
4015     comm = dd->comm;
4016
4017     if (comm->bBondComm)
4018     {
4019         /* Communicate atoms beyond the cut-off for bonded interactions */
4020         comm = dd->comm;
4021
4022         comm->cglink = make_charge_group_links(mtop, dd, cginfo_mb);
4023
4024         comm->bLocalCG = init_bLocalCG(mtop);
4025     }
4026     else
4027     {
4028         /* Only communicate atoms based on cut-off */
4029         comm->cglink   = nullptr;
4030         comm->bLocalCG = nullptr;
4031     }
4032 }
4033
4034 static void print_dd_settings(FILE *fplog, gmx_domdec_t *dd,
4035                               const gmx_mtop_t *mtop, const t_inputrec *ir,
4036                               gmx_bool bDynLoadBal, real dlb_scale,
4037                               const gmx_ddbox_t *ddbox)
4038 {
4039     gmx_domdec_comm_t *comm;
4040     int                d;
4041     ivec               np;
4042     real               limit, shrink;
4043     char               buf[64];
4044
4045     if (fplog == nullptr)
4046     {
4047         return;
4048     }
4049
4050     comm = dd->comm;
4051
4052     if (bDynLoadBal)
4053     {
4054         fprintf(fplog, "The maximum number of communication pulses is:");
4055         for (d = 0; d < dd->ndim; d++)
4056         {
4057             fprintf(fplog, " %c %d", dim2char(dd->dim[d]), comm->cd[d].np_dlb);
4058         }
4059         fprintf(fplog, "\n");
4060         fprintf(fplog, "The minimum size for domain decomposition cells is %.3f nm\n", comm->cellsize_limit);
4061         fprintf(fplog, "The requested allowed shrink of DD cells (option -dds) is: %.2f\n", dlb_scale);
4062         fprintf(fplog, "The allowed shrink of domain decomposition cells is:");
4063         for (d = 0; d < DIM; d++)
4064         {
4065             if (dd->nc[d] > 1)
4066             {
4067                 if (d >= ddbox->npbcdim && dd->nc[d] == 2)
4068                 {
4069                     shrink = 0;
4070                 }
4071                 else
4072                 {
4073                     shrink =
4074                         comm->cellsize_min_dlb[d]/
4075                         (ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
4076                 }
4077                 fprintf(fplog, " %c %.2f", dim2char(d), shrink);
4078             }
4079         }
4080         fprintf(fplog, "\n");
4081     }
4082     else
4083     {
4084         set_dd_cell_sizes_slb(dd, ddbox, setcellsizeslbPULSE_ONLY, np);
4085         fprintf(fplog, "The initial number of communication pulses is:");
4086         for (d = 0; d < dd->ndim; d++)
4087         {
4088             fprintf(fplog, " %c %d", dim2char(dd->dim[d]), np[dd->dim[d]]);
4089         }
4090         fprintf(fplog, "\n");
4091         fprintf(fplog, "The initial domain decomposition cell size is:");
4092         for (d = 0; d < DIM; d++)
4093         {
4094             if (dd->nc[d] > 1)
4095             {
4096                 fprintf(fplog, " %c %.2f nm",
4097                         dim2char(d), dd->comm->cellsize_min[d]);
4098             }
4099         }
4100         fprintf(fplog, "\n\n");
4101     }
4102
4103     gmx_bool bInterCGVsites = count_intercg_vsites(mtop);
4104
4105     if (comm->bInterCGBondeds ||
4106         bInterCGVsites ||
4107         dd->bInterCGcons || dd->bInterCGsettles)
4108     {
4109         fprintf(fplog, "The maximum allowed distance for charge groups involved in interactions is:\n");
4110         fprintf(fplog, "%40s  %-7s %6.3f nm\n",
4111                 "non-bonded interactions", "", comm->cutoff);
4112
4113         if (bDynLoadBal)
4114         {
4115             limit = dd->comm->cellsize_limit;
4116         }
4117         else
4118         {
4119             if (dynamic_dd_box(ddbox, ir))
4120             {
4121                 fprintf(fplog, "(the following are initial values, they could change due to box deformation)\n");
4122             }
4123             limit = dd->comm->cellsize_min[XX];
4124             for (d = 1; d < DIM; d++)
4125             {
4126                 limit = std::min(limit, dd->comm->cellsize_min[d]);
4127             }
4128         }
4129
4130         if (comm->bInterCGBondeds)
4131         {
4132             fprintf(fplog, "%40s  %-7s %6.3f nm\n",
4133                     "two-body bonded interactions", "(-rdd)",
4134                     std::max(comm->cutoff, comm->cutoff_mbody));
4135             fprintf(fplog, "%40s  %-7s %6.3f nm\n",
4136                     "multi-body bonded interactions", "(-rdd)",
4137                     (comm->bBondComm || isDlbOn(dd->comm)) ? comm->cutoff_mbody : std::min(comm->cutoff, limit));
4138         }
4139         if (bInterCGVsites)
4140         {
4141             fprintf(fplog, "%40s  %-7s %6.3f nm\n",
4142                     "virtual site constructions", "(-rcon)", limit);
4143         }
4144         if (dd->bInterCGcons || dd->bInterCGsettles)
4145         {
4146             sprintf(buf, "atoms separated by up to %d constraints",
4147                     1+ir->nProjOrder);
4148             fprintf(fplog, "%40s  %-7s %6.3f nm\n",
4149                     buf, "(-rcon)", limit);
4150         }
4151         fprintf(fplog, "\n");
4152     }
4153
4154     fflush(fplog);
4155 }
4156
4157 static void set_cell_limits_dlb(gmx_domdec_t      *dd,
4158                                 real               dlb_scale,
4159                                 const t_inputrec  *ir,
4160                                 const gmx_ddbox_t *ddbox)
4161 {
4162     gmx_domdec_comm_t *comm;
4163     int                d, dim, npulse, npulse_d_max, npulse_d;
4164     gmx_bool           bNoCutOff;
4165
4166     comm = dd->comm;
4167
4168     bNoCutOff = (ir->rvdw == 0 || ir->rcoulomb == 0);
4169
4170     /* Determine the maximum number of comm. pulses in one dimension */
4171
4172     comm->cellsize_limit = std::max(comm->cellsize_limit, comm->cutoff_mbody);
4173
4174     /* Determine the maximum required number of grid pulses */
4175     if (comm->cellsize_limit >= comm->cutoff)
4176     {
4177         /* Only a single pulse is required */
4178         npulse = 1;
4179     }
4180     else if (!bNoCutOff && comm->cellsize_limit > 0)
4181     {
4182         /* We round down slightly here to avoid overhead due to the latency
4183          * of extra communication calls when the cut-off
4184          * would be only slightly longer than the cell size.
4185          * Later cellsize_limit is redetermined,
4186          * so we can not miss interactions due to this rounding.
4187          */
4188         npulse = (int)(0.96 + comm->cutoff/comm->cellsize_limit);
4189     }
4190     else
4191     {
4192         /* There is no cell size limit */
4193         npulse = std::max(dd->nc[XX]-1, std::max(dd->nc[YY]-1, dd->nc[ZZ]-1));
4194     }
4195
4196     if (!bNoCutOff && npulse > 1)
4197     {
4198         /* See if we can do with less pulses, based on dlb_scale */
4199         npulse_d_max = 0;
4200         for (d = 0; d < dd->ndim; d++)
4201         {
4202             dim      = dd->dim[d];
4203             npulse_d = (int)(1 + dd->nc[dim]*comm->cutoff
4204                              /(ddbox->box_size[dim]*ddbox->skew_fac[dim]*dlb_scale));
4205             npulse_d_max = std::max(npulse_d_max, npulse_d);
4206         }
4207         npulse = std::min(npulse, npulse_d_max);
4208     }
4209
4210     /* This env var can override npulse */
4211     d = dd_getenv(debug, "GMX_DD_NPULSE", 0);
4212     if (d > 0)
4213     {
4214         npulse = d;
4215     }
4216
4217     comm->maxpulse       = 1;
4218     comm->bVacDLBNoLimit = (ir->ePBC == epbcNONE);
4219     for (d = 0; d < dd->ndim; d++)
4220     {
4221         comm->cd[d].np_dlb    = std::min(npulse, dd->nc[dd->dim[d]]-1);
4222         comm->maxpulse        = std::max(comm->maxpulse, comm->cd[d].np_dlb);
4223         if (comm->cd[d].np_dlb < dd->nc[dd->dim[d]]-1)
4224         {
4225             comm->bVacDLBNoLimit = FALSE;
4226         }
4227     }
4228
4229     /* cellsize_limit is set for LINCS in init_domain_decomposition */
4230     if (!comm->bVacDLBNoLimit)
4231     {
4232         comm->cellsize_limit = std::max(comm->cellsize_limit,
4233                                         comm->cutoff/comm->maxpulse);
4234     }
4235     comm->cellsize_limit = std::max(comm->cellsize_limit, comm->cutoff_mbody);
4236     /* Set the minimum cell size for each DD dimension */
4237     for (d = 0; d < dd->ndim; d++)
4238     {
4239         if (comm->bVacDLBNoLimit ||
4240             comm->cd[d].np_dlb*comm->cellsize_limit >= comm->cutoff)
4241         {
4242             comm->cellsize_min_dlb[dd->dim[d]] = comm->cellsize_limit;
4243         }
4244         else
4245         {
4246             comm->cellsize_min_dlb[dd->dim[d]] =
4247                 comm->cutoff/comm->cd[d].np_dlb;
4248         }
4249     }
4250     if (comm->cutoff_mbody <= 0)
4251     {
4252         comm->cutoff_mbody = std::min(comm->cutoff, comm->cellsize_limit);
4253     }
4254     if (isDlbOn(comm))
4255     {
4256         set_dlb_limits(dd);
4257     }
4258 }
4259
4260 gmx_bool dd_bonded_molpbc(const gmx_domdec_t *dd, int ePBC)
4261 {
4262     /* If each molecule is a single charge group
4263      * or we use domain decomposition for each periodic dimension,
4264      * we do not need to take pbc into account for the bonded interactions.
4265      */
4266     return (ePBC != epbcNONE && dd->comm->bInterCGBondeds &&
4267             !(dd->nc[XX] > 1 &&
4268               dd->nc[YY] > 1 &&
4269               (dd->nc[ZZ] > 1 || ePBC == epbcXY)));
4270 }
4271
4272 /*! \brief Sets grid size limits and PP-PME setup, prints settings to log */
4273 static void set_ddgrid_parameters(FILE *fplog, gmx_domdec_t *dd, real dlb_scale,
4274                                   const gmx_mtop_t *mtop, const t_inputrec *ir,
4275                                   const gmx_ddbox_t *ddbox)
4276 {
4277     gmx_domdec_comm_t *comm;
4278     int                natoms_tot;
4279     real               vol_frac;
4280
4281     comm = dd->comm;
4282
4283     if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
4284     {
4285         init_ddpme(dd, &comm->ddpme[0], 0);
4286         if (comm->npmedecompdim >= 2)
4287         {
4288             init_ddpme(dd, &comm->ddpme[1], 1);
4289         }
4290     }
4291     else
4292     {
4293         comm->npmenodes = 0;
4294         if (dd->pme_nodeid >= 0)
4295         {
4296             gmx_fatal_collective(FARGS, dd->mpi_comm_all, DDMASTER(dd),
4297                                  "Can not have separate PME ranks without PME electrostatics");
4298         }
4299     }
4300
4301     if (debug)
4302     {
4303         fprintf(debug, "The DD cut-off is %f\n", comm->cutoff);
4304     }
4305     if (!isDlbDisabled(comm))
4306     {
4307         set_cell_limits_dlb(dd, dlb_scale, ir, ddbox);
4308     }
4309
4310     print_dd_settings(fplog, dd, mtop, ir, isDlbOn(comm), dlb_scale, ddbox);
4311     if (comm->dlbState == edlbsOffCanTurnOn)
4312     {
4313         if (fplog)
4314         {
4315             fprintf(fplog, "When dynamic load balancing gets turned on, these settings will change to:\n");
4316         }
4317         print_dd_settings(fplog, dd, mtop, ir, TRUE, dlb_scale, ddbox);
4318     }
4319
4320     if (ir->ePBC == epbcNONE)
4321     {
4322         vol_frac = 1 - 1/(double)dd->nnodes;
4323     }
4324     else
4325     {
4326         vol_frac =
4327             (1 + comm_box_frac(dd->nc, comm->cutoff, ddbox))/(double)dd->nnodes;
4328     }
4329     if (debug)
4330     {
4331         fprintf(debug, "Volume fraction for all DD zones: %f\n", vol_frac);
4332     }
4333     natoms_tot = comm->cgs_gl.index[comm->cgs_gl.nr];
4334
4335     dd->ga2la = ga2la_init(natoms_tot, static_cast<int>(vol_frac*natoms_tot));
4336 }
4337
4338 /*! \brief Set some important DD parameters that can be modified by env.vars */
4339 static void set_dd_envvar_options(FILE *fplog, gmx_domdec_t *dd, int rank_mysim)
4340 {
4341     gmx_domdec_comm_t *comm = dd->comm;
4342
4343     dd->bSendRecv2      = dd_getenv(fplog, "GMX_DD_USE_SENDRECV2", 0);
4344     comm->dlb_scale_lim = dd_getenv(fplog, "GMX_DLB_MAX_BOX_SCALING", 10);
4345     comm->eFlop         = dd_getenv(fplog, "GMX_DLB_BASED_ON_FLOPS", 0);
4346     int recload         = dd_getenv(fplog, "GMX_DD_RECORD_LOAD", 1);
4347     comm->nstDDDump     = dd_getenv(fplog, "GMX_DD_NST_DUMP", 0);
4348     comm->nstDDDumpGrid = dd_getenv(fplog, "GMX_DD_NST_DUMP_GRID", 0);
4349     comm->DD_debug      = dd_getenv(fplog, "GMX_DD_DEBUG", 0);
4350
4351     if (dd->bSendRecv2 && fplog)
4352     {
4353         fprintf(fplog, "Will use two sequential MPI_Sendrecv calls instead of two simultaneous non-blocking MPI_Irecv and MPI_Isend pairs for constraint and vsite communication\n");
4354     }
4355
4356     if (comm->eFlop)
4357     {
4358         if (fplog)
4359         {
4360             fprintf(fplog, "Will load balance based on FLOP count\n");
4361         }
4362         if (comm->eFlop > 1)
4363         {
4364             srand(1 + rank_mysim);
4365         }
4366         comm->bRecordLoad = TRUE;
4367     }
4368     else
4369     {
4370         comm->bRecordLoad = (wallcycle_have_counter() && recload > 0);
4371     }
4372 }
4373
4374 DomdecOptions::DomdecOptions() :
4375     checkBondedInteractions(TRUE),
4376     useBondedCommunication(TRUE),
4377     numPmeRanks(-1),
4378     rankOrder(DdRankOrder::pp_pme),
4379     minimumCommunicationRange(0),
4380     constraintCommunicationRange(0),
4381     dlbOption(DlbOption::turnOnWhenUseful),
4382     dlbScaling(0.8),
4383     cellSizeX(nullptr),
4384     cellSizeY(nullptr),
4385     cellSizeZ(nullptr)
4386 {
4387     clear_ivec(numCells);
4388 }
4389
4390 gmx_domdec_t *init_domain_decomposition(FILE *fplog, t_commrec *cr,
4391                                         const DomdecOptions &options,
4392                                         const MdrunOptions &mdrunOptions,
4393                                         const gmx_mtop_t *mtop,
4394                                         const t_inputrec *ir,
4395                                         const matrix box,
4396                                         gmx::ArrayRef<const gmx::RVec> xGlobal,
4397                                         gmx::LocalAtomSetManager *atomSets)
4398 {
4399     gmx_domdec_t      *dd;
4400
4401     if (fplog)
4402     {
4403         fprintf(fplog,
4404                 "\nInitializing Domain Decomposition on %d ranks\n", cr->nnodes);
4405     }
4406
4407     dd = new gmx_domdec_t;
4408
4409     dd->comm = init_dd_comm();
4410
4411     /* Initialize DD paritioning counters */
4412     dd->comm->partition_step = INT_MIN;
4413     dd->ddp_count            = 0;
4414
4415     set_dd_envvar_options(fplog, dd, cr->nodeid);
4416
4417     gmx_ddbox_t ddbox = {0};
4418     set_dd_limits_and_grid(fplog, cr, dd, options, mdrunOptions,
4419                            mtop, ir,
4420                            box, xGlobal,
4421                            &ddbox);
4422
4423     make_dd_communicators(fplog, cr, dd, options.rankOrder);
4424
4425     if (thisRankHasDuty(cr, DUTY_PP))
4426     {
4427         set_ddgrid_parameters(fplog, dd, options.dlbScaling, mtop, ir, &ddbox);
4428
4429         setup_neighbor_relations(dd);
4430     }
4431
4432     /* Set overallocation to avoid frequent reallocation of arrays */
4433     set_over_alloc_dd(TRUE);
4434
4435     clear_dd_cycle_counts(dd);
4436
4437     dd->atomSets = atomSets;
4438
4439     return dd;
4440 }
4441
4442 static gmx_bool test_dd_cutoff(t_commrec *cr,
4443                                t_state *state, const t_inputrec *ir,
4444                                real cutoff_req)
4445 {
4446     gmx_domdec_t *dd;
4447     gmx_ddbox_t   ddbox;
4448     int           d, dim, np;
4449     real          inv_cell_size;
4450     int           LocallyLimited;
4451
4452     dd = cr->dd;
4453
4454     set_ddbox(dd, false, ir, state->box, true, state->x, &ddbox);
4455
4456     LocallyLimited = 0;
4457
4458     for (d = 0; d < dd->ndim; d++)
4459     {
4460         dim = dd->dim[d];
4461
4462         inv_cell_size = DD_CELL_MARGIN*dd->nc[dim]/ddbox.box_size[dim];
4463         if (dynamic_dd_box(&ddbox, ir))
4464         {
4465             inv_cell_size *= DD_PRES_SCALE_MARGIN;
4466         }
4467
4468         np = 1 + (int)(cutoff_req*inv_cell_size*ddbox.skew_fac[dim]);
4469
4470         if (!isDlbDisabled(dd->comm) && (dim < ddbox.npbcdim) && (dd->comm->cd[d].np_dlb > 0))
4471         {
4472             if (np > dd->comm->cd[d].np_dlb)
4473             {
4474                 return FALSE;
4475             }
4476
4477             /* If a current local cell size is smaller than the requested
4478              * cut-off, we could still fix it, but this gets very complicated.
4479              * Without fixing here, we might actually need more checks.
4480              */
4481             if ((dd->comm->cell_x1[dim] - dd->comm->cell_x0[dim])*ddbox.skew_fac[dim]*dd->comm->cd[d].np_dlb < cutoff_req)
4482             {
4483                 LocallyLimited = 1;
4484             }
4485         }
4486     }
4487
4488     if (!isDlbDisabled(dd->comm))
4489     {
4490         /* If DLB is not active yet, we don't need to check the grid jumps.
4491          * Actually we shouldn't, because then the grid jump data is not set.
4492          */
4493         if (isDlbOn(dd->comm) &&
4494             check_grid_jump(0, dd, cutoff_req, &ddbox, FALSE))
4495         {
4496             LocallyLimited = 1;
4497         }
4498
4499         gmx_sumi(1, &LocallyLimited, cr);
4500
4501         if (LocallyLimited > 0)
4502         {
4503             return FALSE;
4504         }
4505     }
4506
4507     return TRUE;
4508 }
4509
4510 gmx_bool change_dd_cutoff(t_commrec *cr, t_state *state, const t_inputrec *ir,
4511                           real cutoff_req)
4512 {
4513     gmx_bool bCutoffAllowed;
4514
4515     bCutoffAllowed = test_dd_cutoff(cr, state, ir, cutoff_req);
4516
4517     if (bCutoffAllowed)
4518     {
4519         cr->dd->comm->cutoff = cutoff_req;
4520     }
4521
4522     return bCutoffAllowed;
4523 }
4524
4525 void set_dd_dlb_max_cutoff(t_commrec *cr, real cutoff)
4526 {
4527     gmx_domdec_comm_t *comm;
4528
4529     comm = cr->dd->comm;
4530
4531     /* Turn on the DLB limiting (might have been on already) */
4532     comm->bPMELoadBalDLBLimits = TRUE;
4533
4534     /* Change the cut-off limit */
4535     comm->PMELoadBal_max_cutoff = cutoff;
4536
4537     if (debug)
4538     {
4539         fprintf(debug, "PME load balancing set a limit to the DLB staggering such that a %f cut-off will continue to fit\n",
4540                 comm->PMELoadBal_max_cutoff);
4541     }
4542 }
4543
4544 /* Sets whether we should later check the load imbalance data, so that
4545  * we can trigger dynamic load balancing if enough imbalance has
4546  * arisen.
4547  *
4548  * Used after PME load balancing unlocks DLB, so that the check
4549  * whether DLB will be useful can happen immediately.
4550  */
4551 static void dd_dlb_set_should_check_whether_to_turn_dlb_on(gmx_domdec_t *dd, gmx_bool bValue)
4552 {
4553     if (dd->comm->dlbState == edlbsOffCanTurnOn)
4554     {
4555         dd->comm->bCheckWhetherToTurnDlbOn = bValue;
4556
4557         if (bValue == TRUE)
4558         {
4559             /* Store the DD partitioning count, so we can ignore cycle counts
4560              * over the next nstlist steps, which are often slower.
4561              */
4562             dd->comm->ddPartioningCountFirstDlbOff = dd->ddp_count;
4563         }
4564     }
4565 }
4566
4567 /* Returns if we should check whether there has been enough load
4568  * imbalance to trigger dynamic load balancing.
4569  */
4570 static gmx_bool dd_dlb_get_should_check_whether_to_turn_dlb_on(gmx_domdec_t *dd)
4571 {
4572     if (dd->comm->dlbState != edlbsOffCanTurnOn)
4573     {
4574         return FALSE;
4575     }
4576
4577     if (dd->ddp_count <= dd->comm->ddPartioningCountFirstDlbOff)
4578     {
4579         /* We ignore the first nstlist steps at the start of the run
4580          * or after PME load balancing or after turning DLB off, since
4581          * these often have extra allocation or cache miss overhead.
4582          */
4583         return FALSE;
4584     }
4585
4586     if (dd->comm->cycl_n[ddCyclStep] == 0)
4587     {
4588         /* We can have zero timed steps when dd_partition_system is called
4589          * more than once at the same step, e.g. with replica exchange.
4590          * Turning on DLB would trigger an assertion failure later, but is
4591          * also useless right after exchanging replicas.
4592          */
4593         return FALSE;
4594     }
4595
4596     /* We should check whether we should use DLB directly after
4597      * unlocking DLB. */
4598     if (dd->comm->bCheckWhetherToTurnDlbOn)
4599     {
4600         /* This flag was set when the PME load-balancing routines
4601            unlocked DLB, and should now be cleared. */
4602         dd_dlb_set_should_check_whether_to_turn_dlb_on(dd, FALSE);
4603         return TRUE;
4604     }
4605     /* We check whether we should use DLB every c_checkTurnDlbOnInterval
4606      * partitionings (we do not do this every partioning, so that we
4607      * avoid excessive communication). */
4608     return dd->comm->n_load_have % c_checkTurnDlbOnInterval == c_checkTurnDlbOnInterval - 1;
4609 }
4610
4611 gmx_bool dd_dlb_is_on(const gmx_domdec_t *dd)
4612 {
4613     return isDlbOn(dd->comm);
4614 }
4615
4616 gmx_bool dd_dlb_is_locked(const gmx_domdec_t *dd)
4617 {
4618     return (dd->comm->dlbState == edlbsOffTemporarilyLocked);
4619 }
4620
4621 void dd_dlb_lock(gmx_domdec_t *dd)
4622 {
4623     /* We can only lock the DLB when it is set to auto, otherwise don't do anything */
4624     if (dd->comm->dlbState == edlbsOffCanTurnOn)
4625     {
4626         dd->comm->dlbState = edlbsOffTemporarilyLocked;
4627     }
4628 }
4629
4630 void dd_dlb_unlock(gmx_domdec_t *dd)
4631 {
4632     /* We can only lock the DLB when it is set to auto, otherwise don't do anything */
4633     if (dd->comm->dlbState == edlbsOffTemporarilyLocked)
4634     {
4635         dd->comm->dlbState = edlbsOffCanTurnOn;
4636         dd_dlb_set_should_check_whether_to_turn_dlb_on(dd, TRUE);
4637     }
4638 }
4639
4640 static void merge_cg_buffers(int ncell,
4641                              gmx_domdec_comm_dim_t *cd, int pulse,
4642                              int  *ncg_cell,
4643                              gmx::ArrayRef<int> index_gl,
4644                              const int  *recv_i,
4645                              rvec *cg_cm,    rvec *recv_vr,
4646                              gmx::ArrayRef<int> cgindex,
4647                              cginfo_mb_t *cginfo_mb, int *cginfo)
4648 {
4649     gmx_domdec_ind_t *ind, *ind_p;
4650     int               p, cell, c, cg, cg0, cg1, cg_gl, nat;
4651     int               shift, shift_at;
4652
4653     ind = &cd->ind[pulse];
4654
4655     /* First correct the already stored data */
4656     shift = ind->nrecv[ncell];
4657     for (cell = ncell-1; cell >= 0; cell--)
4658     {
4659         shift -= ind->nrecv[cell];
4660         if (shift > 0)
4661         {
4662             /* Move the cg's present from previous grid pulses */
4663             cg0                = ncg_cell[ncell+cell];
4664             cg1                = ncg_cell[ncell+cell+1];
4665             cgindex[cg1+shift] = cgindex[cg1];
4666             for (cg = cg1-1; cg >= cg0; cg--)
4667             {
4668                 index_gl[cg+shift] = index_gl[cg];
4669                 copy_rvec(cg_cm[cg], cg_cm[cg+shift]);
4670                 cgindex[cg+shift] = cgindex[cg];
4671                 cginfo[cg+shift]  = cginfo[cg];
4672             }
4673             /* Correct the already stored send indices for the shift */
4674             for (p = 1; p <= pulse; p++)
4675             {
4676                 ind_p = &cd->ind[p];
4677                 cg0   = 0;
4678                 for (c = 0; c < cell; c++)
4679                 {
4680                     cg0 += ind_p->nsend[c];
4681                 }
4682                 cg1 = cg0 + ind_p->nsend[cell];
4683                 for (cg = cg0; cg < cg1; cg++)
4684                 {
4685                     ind_p->index[cg] += shift;
4686                 }
4687             }
4688         }
4689     }
4690
4691     /* Merge in the communicated buffers */
4692     shift    = 0;
4693     shift_at = 0;
4694     cg0      = 0;
4695     for (cell = 0; cell < ncell; cell++)
4696     {
4697         cg1 = ncg_cell[ncell+cell+1] + shift;
4698         if (shift_at > 0)
4699         {
4700             /* Correct the old cg indices */
4701             for (cg = ncg_cell[ncell+cell]; cg < cg1; cg++)
4702             {
4703                 cgindex[cg+1] += shift_at;
4704             }
4705         }
4706         for (cg = 0; cg < ind->nrecv[cell]; cg++)
4707         {
4708             /* Copy this charge group from the buffer */
4709             index_gl[cg1] = recv_i[cg0];
4710             copy_rvec(recv_vr[cg0], cg_cm[cg1]);
4711             /* Add it to the cgindex */
4712             cg_gl          = index_gl[cg1];
4713             cginfo[cg1]    = ddcginfo(cginfo_mb, cg_gl);
4714             nat            = GET_CGINFO_NATOMS(cginfo[cg1]);
4715             cgindex[cg1+1] = cgindex[cg1] + nat;
4716             cg0++;
4717             cg1++;
4718             shift_at += nat;
4719         }
4720         shift                 += ind->nrecv[cell];
4721         ncg_cell[ncell+cell+1] = cg1;
4722     }
4723 }
4724
4725 static void make_cell2at_index(gmx_domdec_comm_dim_t        *cd,
4726                                int                           nzone,
4727                                int                           atomGroupStart,
4728                                const gmx::RangePartitioning &atomGroups)
4729 {
4730     /* Store the atom block boundaries for easy copying of communication buffers
4731      */
4732     int g = atomGroupStart;
4733     for (int zone = 0; zone < nzone; zone++)
4734     {
4735         for (gmx_domdec_ind_t &ind : cd->ind)
4736         {
4737             const auto range    = atomGroups.subRange(g, g + ind.nrecv[zone]);
4738             ind.cell2at0[zone]  = range.begin();
4739             ind.cell2at1[zone]  = range.end();
4740             g                  += ind.nrecv[zone];
4741         }
4742     }
4743 }
4744
4745 static gmx_bool missing_link(t_blocka *link, int cg_gl, const char *bLocalCG)
4746 {
4747     int      i;
4748     gmx_bool bMiss;
4749
4750     bMiss = FALSE;
4751     for (i = link->index[cg_gl]; i < link->index[cg_gl+1]; i++)
4752     {
4753         if (!bLocalCG[link->a[i]])
4754         {
4755             bMiss = TRUE;
4756         }
4757     }
4758
4759     return bMiss;
4760 }
4761
4762 /* Domain corners for communication, a maximum of 4 i-zones see a j domain */
4763 typedef struct {
4764     real c[DIM][4]; /* the corners for the non-bonded communication */
4765     real cr0;       /* corner for rounding */
4766     real cr1[4];    /* corners for rounding */
4767     real bc[DIM];   /* corners for bounded communication */
4768     real bcr1;      /* corner for rounding for bonded communication */
4769 } dd_corners_t;
4770
4771 /* Determine the corners of the domain(s) we are communicating with */
4772 static void
4773 set_dd_corners(const gmx_domdec_t *dd,
4774                int dim0, int dim1, int dim2,
4775                gmx_bool bDistMB,
4776                dd_corners_t *c)
4777 {
4778     const gmx_domdec_comm_t  *comm;
4779     const gmx_domdec_zones_t *zones;
4780     int i, j;
4781
4782     comm = dd->comm;
4783
4784     zones = &comm->zones;
4785
4786     /* Keep the compiler happy */
4787     c->cr0  = 0;
4788     c->bcr1 = 0;
4789
4790     /* The first dimension is equal for all cells */
4791     c->c[0][0] = comm->cell_x0[dim0];
4792     if (bDistMB)
4793     {
4794         c->bc[0] = c->c[0][0];
4795     }
4796     if (dd->ndim >= 2)
4797     {
4798         dim1 = dd->dim[1];
4799         /* This cell row is only seen from the first row */
4800         c->c[1][0] = comm->cell_x0[dim1];
4801         /* All rows can see this row */
4802         c->c[1][1] = comm->cell_x0[dim1];
4803         if (isDlbOn(dd->comm))
4804         {
4805             c->c[1][1] = std::max(comm->cell_x0[dim1], comm->zone_d1[1].mch0);
4806             if (bDistMB)
4807             {
4808                 /* For the multi-body distance we need the maximum */
4809                 c->bc[1] = std::max(comm->cell_x0[dim1], comm->zone_d1[1].p1_0);
4810             }
4811         }
4812         /* Set the upper-right corner for rounding */
4813         c->cr0 = comm->cell_x1[dim0];
4814
4815         if (dd->ndim >= 3)
4816         {
4817             dim2 = dd->dim[2];
4818             for (j = 0; j < 4; j++)
4819             {
4820                 c->c[2][j] = comm->cell_x0[dim2];
4821             }
4822             if (isDlbOn(dd->comm))
4823             {
4824                 /* Use the maximum of the i-cells that see a j-cell */
4825                 for (i = 0; i < zones->nizone; i++)
4826                 {
4827                     for (j = zones->izone[i].j0; j < zones->izone[i].j1; j++)
4828                     {
4829                         if (j >= 4)
4830                         {
4831                             c->c[2][j-4] =
4832                                 std::max(c->c[2][j-4],
4833                                          comm->zone_d2[zones->shift[i][dim0]][zones->shift[i][dim1]].mch0);
4834                         }
4835                     }
4836                 }
4837                 if (bDistMB)
4838                 {
4839                     /* For the multi-body distance we need the maximum */
4840                     c->bc[2] = comm->cell_x0[dim2];
4841                     for (i = 0; i < 2; i++)
4842                     {
4843                         for (j = 0; j < 2; j++)
4844                         {
4845                             c->bc[2] = std::max(c->bc[2], comm->zone_d2[i][j].p1_0);
4846                         }
4847                     }
4848                 }
4849             }
4850
4851             /* Set the upper-right corner for rounding */
4852             /* Cell (0,0,0) and cell (1,0,0) can see cell 4 (0,1,1)
4853              * Only cell (0,0,0) can see cell 7 (1,1,1)
4854              */
4855             c->cr1[0] = comm->cell_x1[dim1];
4856             c->cr1[3] = comm->cell_x1[dim1];
4857             if (isDlbOn(dd->comm))
4858             {
4859                 c->cr1[0] = std::max(comm->cell_x1[dim1], comm->zone_d1[1].mch1);
4860                 if (bDistMB)
4861                 {
4862                     /* For the multi-body distance we need the maximum */
4863                     c->bcr1 = std::max(comm->cell_x1[dim1], comm->zone_d1[1].p1_1);
4864                 }
4865             }
4866         }
4867     }
4868 }
4869
4870 /* Add the atom groups we need to send in this pulse from this zone to
4871  * \p localAtomGroups and \p work
4872  */
4873 static void
4874 get_zone_pulse_cgs(gmx_domdec_t *dd,
4875                    int zonei, int zone,
4876                    int cg0, int cg1,
4877                    gmx::ArrayRef<const int> globalAtomGroupIndices,
4878                    const gmx::RangePartitioning &atomGroups,
4879                    int dim, int dim_ind,
4880                    int dim0, int dim1, int dim2,
4881                    real r_comm2, real r_bcomm2,
4882                    matrix box,
4883                    bool distanceIsTriclinic,
4884                    rvec *normal,
4885                    real skew_fac2_d, real skew_fac_01,
4886                    rvec *v_d, rvec *v_0, rvec *v_1,
4887                    const dd_corners_t *c,
4888                    const rvec sf2_round,
4889                    gmx_bool bDistBonded,
4890                    gmx_bool bBondComm,
4891                    gmx_bool bDist2B,
4892                    gmx_bool bDistMB,
4893                    rvec *cg_cm,
4894                    const int *cginfo,
4895                    std::vector<int>     *localAtomGroups,
4896                    dd_comm_setup_work_t *work)
4897 {
4898     gmx_domdec_comm_t *comm;
4899     gmx_bool           bScrew;
4900     gmx_bool           bDistMB_pulse;
4901     int                cg, i;
4902     real               r2, rb2, r, tric_sh;
4903     rvec               rn, rb;
4904     int                dimd;
4905     int                nsend_z, nat;
4906
4907     comm = dd->comm;
4908
4909     bScrew = (dd->bScrewPBC && dim == XX);
4910
4911     bDistMB_pulse = (bDistMB && bDistBonded);
4912
4913     /* Unpack the work data */
4914     std::vector<int>       &ibuf = work->atomGroupBuffer;
4915     std::vector<gmx::RVec> &vbuf = work->positionBuffer;
4916     nsend_z                      = 0;
4917     nat                          = work->nat;
4918
4919     for (cg = cg0; cg < cg1; cg++)
4920     {
4921         r2  = 0;
4922         rb2 = 0;
4923         if (!distanceIsTriclinic)
4924         {
4925             /* Rectangular direction, easy */
4926             r = cg_cm[cg][dim] - c->c[dim_ind][zone];
4927             if (r > 0)
4928             {
4929                 r2 += r*r;
4930             }
4931             if (bDistMB_pulse)
4932             {
4933                 r = cg_cm[cg][dim] - c->bc[dim_ind];
4934                 if (r > 0)
4935                 {
4936                     rb2 += r*r;
4937                 }
4938             }
4939             /* Rounding gives at most a 16% reduction
4940              * in communicated atoms
4941              */
4942             if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
4943             {
4944                 r = cg_cm[cg][dim0] - c->cr0;
4945                 /* This is the first dimension, so always r >= 0 */
4946                 r2 += r*r;
4947                 if (bDistMB_pulse)
4948                 {
4949                     rb2 += r*r;
4950                 }
4951             }
4952             if (dim_ind == 2 && (zonei == 2 || zonei == 3))
4953             {
4954                 r = cg_cm[cg][dim1] - c->cr1[zone];
4955                 if (r > 0)
4956                 {
4957                     r2 += r*r;
4958                 }
4959                 if (bDistMB_pulse)
4960                 {
4961                     r = cg_cm[cg][dim1] - c->bcr1;
4962                     if (r > 0)
4963                     {
4964                         rb2 += r*r;
4965                     }
4966                 }
4967             }
4968         }
4969         else
4970         {
4971             /* Triclinic direction, more complicated */
4972             clear_rvec(rn);
4973             clear_rvec(rb);
4974             /* Rounding, conservative as the skew_fac multiplication
4975              * will slightly underestimate the distance.
4976              */
4977             if (dim_ind >= 1 && (zonei == 1 || zonei == 2))
4978             {
4979                 rn[dim0] = cg_cm[cg][dim0] - c->cr0;
4980                 for (i = dim0+1; i < DIM; i++)
4981                 {
4982                     rn[dim0] -= cg_cm[cg][i]*v_0[i][dim0];
4983                 }
4984                 r2 = rn[dim0]*rn[dim0]*sf2_round[dim0];
4985                 if (bDistMB_pulse)
4986                 {
4987                     rb[dim0] = rn[dim0];
4988                     rb2      = r2;
4989                 }
4990                 /* Take care that the cell planes along dim0 might not
4991                  * be orthogonal to those along dim1 and dim2.
4992                  */
4993                 for (i = 1; i <= dim_ind; i++)
4994                 {
4995                     dimd = dd->dim[i];
4996                     if (normal[dim0][dimd] > 0)
4997                     {
4998                         rn[dimd] -= rn[dim0]*normal[dim0][dimd];
4999                         if (bDistMB_pulse)
5000                         {
5001                             rb[dimd] -= rb[dim0]*normal[dim0][dimd];
5002                         }
5003                     }
5004                 }
5005             }
5006             if (dim_ind == 2 && (zonei == 2 || zonei == 3))
5007             {
5008                 rn[dim1] += cg_cm[cg][dim1] - c->cr1[zone];
5009                 tric_sh   = 0;
5010                 for (i = dim1+1; i < DIM; i++)
5011                 {
5012                     tric_sh -= cg_cm[cg][i]*v_1[i][dim1];
5013                 }
5014                 rn[dim1] += tric_sh;
5015                 if (rn[dim1] > 0)
5016                 {
5017                     r2 += rn[dim1]*rn[dim1]*sf2_round[dim1];
5018                     /* Take care of coupling of the distances
5019                      * to the planes along dim0 and dim1 through dim2.
5020                      */
5021                     r2 -= rn[dim0]*rn[dim1]*skew_fac_01;
5022                     /* Take care that the cell planes along dim1
5023                      * might not be orthogonal to that along dim2.
5024                      */
5025                     if (normal[dim1][dim2] > 0)
5026                     {
5027                         rn[dim2] -= rn[dim1]*normal[dim1][dim2];
5028                     }
5029                 }
5030                 if (bDistMB_pulse)
5031                 {
5032                     rb[dim1] +=
5033                         cg_cm[cg][dim1] - c->bcr1 + tric_sh;
5034                     if (rb[dim1] > 0)
5035                     {
5036                         rb2 += rb[dim1]*rb[dim1]*sf2_round[dim1];
5037                         /* Take care of coupling of the distances
5038                          * to the planes along dim0 and dim1 through dim2.
5039                          */
5040                         rb2 -= rb[dim0]*rb[dim1]*skew_fac_01;
5041                         /* Take care that the cell planes along dim1
5042                          * might not be orthogonal to that along dim2.
5043                          */
5044                         if (normal[dim1][dim2] > 0)
5045                         {
5046                             rb[dim2] -= rb[dim1]*normal[dim1][dim2];
5047                         }
5048                     }
5049                 }
5050             }
5051             /* The distance along the communication direction */
5052             rn[dim] += cg_cm[cg][dim] - c->c[dim_ind][zone];
5053             tric_sh  = 0;
5054             for (i = dim+1; i < DIM; i++)
5055             {
5056                 tric_sh -= cg_cm[cg][i]*v_d[i][dim];
5057             }
5058             rn[dim] += tric_sh;
5059             if (rn[dim] > 0)
5060             {
5061                 r2 += rn[dim]*rn[dim]*skew_fac2_d;
5062                 /* Take care of coupling of the distances
5063                  * to the planes along dim0 and dim1 through dim2.
5064                  */
5065                 if (dim_ind == 1 && zonei == 1)
5066                 {
5067                     r2 -= rn[dim0]*rn[dim]*skew_fac_01;
5068                 }
5069             }
5070             if (bDistMB_pulse)
5071             {
5072                 clear_rvec(rb);
5073                 rb[dim] += cg_cm[cg][dim] - c->bc[dim_ind] + tric_sh;
5074                 if (rb[dim] > 0)
5075                 {
5076                     rb2 += rb[dim]*rb[dim]*skew_fac2_d;
5077                     /* Take care of coupling of the distances
5078                      * to the planes along dim0 and dim1 through dim2.
5079                      */
5080                     if (dim_ind == 1 && zonei == 1)
5081                     {
5082                         rb2 -= rb[dim0]*rb[dim]*skew_fac_01;
5083                     }
5084                 }
5085             }
5086         }
5087
5088         if (r2 < r_comm2 ||
5089             (bDistBonded &&
5090              ((bDistMB && rb2 < r_bcomm2) ||
5091               (bDist2B && r2  < r_bcomm2)) &&
5092              (!bBondComm ||
5093               (GET_CGINFO_BOND_INTER(cginfo[cg]) &&
5094                missing_link(comm->cglink, globalAtomGroupIndices[cg],
5095                             comm->bLocalCG)))))
5096         {
5097             /* Store the local and global atom group indices and position */
5098             localAtomGroups->push_back(cg);
5099             ibuf.push_back(globalAtomGroupIndices[cg]);
5100             nsend_z++;
5101
5102             rvec posPbc;
5103             if (dd->ci[dim] == 0)
5104             {
5105                 /* Correct cg_cm for pbc */
5106                 rvec_add(cg_cm[cg], box[dim], posPbc);
5107                 if (bScrew)
5108                 {
5109                     posPbc[YY] = box[YY][YY] - posPbc[YY];
5110                     posPbc[ZZ] = box[ZZ][ZZ] - posPbc[ZZ];
5111                 }
5112             }
5113             else
5114             {
5115                 copy_rvec(cg_cm[cg], posPbc);
5116             }
5117             vbuf.emplace_back(posPbc[XX], posPbc[YY], posPbc[ZZ]);
5118
5119             nat += atomGroups.block(cg).size();
5120         }
5121     }
5122
5123     work->nat        = nat;
5124     work->nsend_zone = nsend_z;
5125 }
5126
5127 static void clearCommSetupData(dd_comm_setup_work_t *work)
5128 {
5129     work->localAtomGroupBuffer.clear();
5130     work->atomGroupBuffer.clear();
5131     work->positionBuffer.clear();
5132     work->nat        = 0;
5133     work->nsend_zone = 0;
5134 }
5135
5136 static void setup_dd_communication(gmx_domdec_t *dd,
5137                                    matrix box, gmx_ddbox_t *ddbox,
5138                                    t_forcerec *fr,
5139                                    t_state *state, PaddedRVecVector *f)
5140 {
5141     int                    dim_ind, dim, dim0, dim1, dim2, dimd, nat_tot;
5142     int                    nzone, nzone_send, zone, zonei, cg0, cg1;
5143     int                    c, i, cg, cg_gl, nrcg;
5144     int                   *zone_cg_range, pos_cg;
5145     gmx_domdec_comm_t     *comm;
5146     gmx_domdec_zones_t    *zones;
5147     gmx_domdec_comm_dim_t *cd;
5148     cginfo_mb_t           *cginfo_mb;
5149     gmx_bool               bBondComm, bDist2B, bDistMB, bDistBonded;
5150     real                   r_comm2, r_bcomm2;
5151     dd_corners_t           corners;
5152     rvec                  *cg_cm, *normal, *v_d, *v_0 = nullptr, *v_1 = nullptr;
5153     real                   skew_fac2_d, skew_fac_01;
5154     rvec                   sf2_round;
5155
5156     if (debug)
5157     {
5158         fprintf(debug, "Setting up DD communication\n");
5159     }
5160
5161     comm  = dd->comm;
5162
5163     if (comm->dth.empty())
5164     {
5165         /* Initialize the thread data.
5166          * This can not be done in init_domain_decomposition,
5167          * as the numbers of threads is determined later.
5168          */
5169         int numThreads = gmx_omp_nthreads_get(emntDomdec);
5170         comm->dth.resize(numThreads);
5171     }
5172
5173     switch (fr->cutoff_scheme)
5174     {
5175         case ecutsGROUP:
5176             cg_cm = fr->cg_cm;
5177             break;
5178         case ecutsVERLET:
5179             cg_cm = as_rvec_array(state->x.data());
5180             break;
5181         default:
5182             gmx_incons("unimplemented");
5183     }
5184
5185     bBondComm = comm->bBondComm;
5186
5187     /* Do we need to determine extra distances for multi-body bondeds? */
5188     bDistMB = (comm->bInterCGMultiBody && isDlbOn(dd->comm) && dd->ndim > 1);
5189
5190     /* Do we need to determine extra distances for only two-body bondeds? */
5191     bDist2B = (bBondComm && !bDistMB);
5192
5193     r_comm2  = gmx::square(comm->cutoff);
5194     r_bcomm2 = gmx::square(comm->cutoff_mbody);
5195
5196     if (debug)
5197     {
5198         fprintf(debug, "bBondComm %s, r_bc %f\n", gmx::boolToString(bBondComm), std::sqrt(r_bcomm2));
5199     }
5200
5201     zones = &comm->zones;
5202
5203     dim0 = dd->dim[0];
5204     dim1 = (dd->ndim >= 2 ? dd->dim[1] : -1);
5205     dim2 = (dd->ndim >= 3 ? dd->dim[2] : -1);
5206
5207     set_dd_corners(dd, dim0, dim1, dim2, bDistMB, &corners);
5208
5209     /* Triclinic stuff */
5210     normal      = ddbox->normal;
5211     skew_fac_01 = 0;
5212     if (dd->ndim >= 2)
5213     {
5214         v_0 = ddbox->v[dim0];
5215         if (ddbox->tric_dir[dim0] && ddbox->tric_dir[dim1])
5216         {
5217             /* Determine the coupling coefficient for the distances
5218              * to the cell planes along dim0 and dim1 through dim2.
5219              * This is required for correct rounding.
5220              */
5221             skew_fac_01 =
5222                 ddbox->v[dim0][dim1+1][dim0]*ddbox->v[dim1][dim1+1][dim1];
5223             if (debug)
5224             {
5225                 fprintf(debug, "\nskew_fac_01 %f\n", skew_fac_01);
5226             }
5227         }
5228     }
5229     if (dd->ndim >= 3)
5230     {
5231         v_1 = ddbox->v[dim1];
5232     }
5233
5234     zone_cg_range = zones->cg_range;
5235     cginfo_mb     = fr->cginfo_mb;
5236
5237     zone_cg_range[0]   = 0;
5238     zone_cg_range[1]   = dd->ncg_home;
5239     comm->zone_ncg1[0] = dd->ncg_home;
5240     pos_cg             = dd->ncg_home;
5241
5242     nat_tot = comm->atomRanges.numHomeAtoms();
5243     nzone   = 1;
5244     for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
5245     {
5246         dim = dd->dim[dim_ind];
5247         cd  = &comm->cd[dim_ind];
5248
5249         /* Check if we need to compute triclinic distances along this dim */
5250         bool distanceIsTriclinic = false;
5251         for (i = 0; i <= dim_ind; i++)
5252         {
5253             if (ddbox->tric_dir[dd->dim[i]])
5254             {
5255                 distanceIsTriclinic = true;
5256             }
5257         }
5258
5259         if (dim >= ddbox->npbcdim && dd->ci[dim] == 0)
5260         {
5261             /* No pbc in this dimension, the first node should not comm. */
5262             nzone_send = 0;
5263         }
5264         else
5265         {
5266             nzone_send = nzone;
5267         }
5268
5269         v_d                = ddbox->v[dim];
5270         skew_fac2_d        = gmx::square(ddbox->skew_fac[dim]);
5271
5272         cd->receiveInPlace = true;
5273         for (int p = 0; p < cd->numPulses(); p++)
5274         {
5275             /* Only atoms communicated in the first pulse are used
5276              * for multi-body bonded interactions or for bBondComm.
5277              */
5278             bDistBonded = ((bDistMB || bDist2B) && p == 0);
5279
5280             gmx_domdec_ind_t *ind = &cd->ind[p];
5281
5282             /* Thread 0 writes in the global index array */
5283             ind->index.clear();
5284             clearCommSetupData(&comm->dth[0]);
5285
5286             for (zone = 0; zone < nzone_send; zone++)
5287             {
5288                 if (dim_ind > 0 && distanceIsTriclinic)
5289                 {
5290                     /* Determine slightly more optimized skew_fac's
5291                      * for rounding.
5292                      * This reduces the number of communicated atoms
5293                      * by about 10% for 3D DD of rhombic dodecahedra.
5294                      */
5295                     for (dimd = 0; dimd < dim; dimd++)
5296                     {
5297                         sf2_round[dimd] = 1;
5298                         if (ddbox->tric_dir[dimd])
5299                         {
5300                             for (i = dd->dim[dimd]+1; i < DIM; i++)
5301                             {
5302                                 /* If we are shifted in dimension i
5303                                  * and the cell plane is tilted forward
5304                                  * in dimension i, skip this coupling.
5305                                  */
5306                                 if (!(zones->shift[nzone+zone][i] &&
5307                                       ddbox->v[dimd][i][dimd] >= 0))
5308                                 {
5309                                     sf2_round[dimd] +=
5310                                         gmx::square(ddbox->v[dimd][i][dimd]);
5311                                 }
5312                             }
5313                             sf2_round[dimd] = 1/sf2_round[dimd];
5314                         }
5315                     }
5316                 }
5317
5318                 zonei = zone_perm[dim_ind][zone];
5319                 if (p == 0)
5320                 {
5321                     /* Here we permutate the zones to obtain a convenient order
5322                      * for neighbor searching
5323                      */
5324                     cg0 = zone_cg_range[zonei];
5325                     cg1 = zone_cg_range[zonei+1];
5326                 }
5327                 else
5328                 {
5329                     /* Look only at the cg's received in the previous grid pulse
5330                      */
5331                     cg1 = zone_cg_range[nzone+zone+1];
5332                     cg0 = cg1 - cd->ind[p-1].nrecv[zone];
5333                 }
5334
5335                 const int numThreads = static_cast<int>(comm->dth.size());
5336 #pragma omp parallel for num_threads(numThreads) schedule(static)
5337                 for (int th = 0; th < numThreads; th++)
5338                 {
5339                     try
5340                     {
5341                         dd_comm_setup_work_t &work = comm->dth[th];
5342
5343                         /* Retain data accumulated into buffers of thread 0 */
5344                         if (th > 0)
5345                         {
5346                             clearCommSetupData(&work);
5347                         }
5348
5349                         int cg0_th = cg0 + ((cg1 - cg0)* th   )/numThreads;
5350                         int cg1_th = cg0 + ((cg1 - cg0)*(th+1))/numThreads;
5351
5352                         /* Get the cg's for this pulse in this zone */
5353                         get_zone_pulse_cgs(dd, zonei, zone, cg0_th, cg1_th,
5354                                            dd->globalAtomGroupIndices,
5355                                            dd->atomGrouping(),
5356                                            dim, dim_ind, dim0, dim1, dim2,
5357                                            r_comm2, r_bcomm2,
5358                                            box, distanceIsTriclinic,
5359                                            normal, skew_fac2_d, skew_fac_01,
5360                                            v_d, v_0, v_1, &corners, sf2_round,
5361                                            bDistBonded, bBondComm,
5362                                            bDist2B, bDistMB,
5363                                            cg_cm, fr->cginfo,
5364                                            th == 0 ? &ind->index : &work.localAtomGroupBuffer,
5365                                            &work);
5366                     }
5367                     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
5368                 } // END
5369
5370                 std::vector<int>       &atomGroups = comm->dth[0].atomGroupBuffer;
5371                 std::vector<gmx::RVec> &positions  = comm->dth[0].positionBuffer;
5372                 ind->nsend[zone]  = comm->dth[0].nsend_zone;
5373                 /* Append data of threads>=1 to the communication buffers */
5374                 for (int th = 1; th < numThreads; th++)
5375                 {
5376                     const dd_comm_setup_work_t &dth = comm->dth[th];
5377
5378                     ind->index.insert(ind->index.end(), dth.localAtomGroupBuffer.begin(), dth.localAtomGroupBuffer.end());
5379                     atomGroups.insert(atomGroups.end(), dth.atomGroupBuffer.begin(), dth.atomGroupBuffer.end());
5380                     positions.insert(positions.end(), dth.positionBuffer.begin(), dth.positionBuffer.end());
5381                     comm->dth[0].nat += dth.nat;
5382                     ind->nsend[zone] += dth.nsend_zone;
5383                 }
5384             }
5385             /* Clear the counts in case we do not have pbc */
5386             for (zone = nzone_send; zone < nzone; zone++)
5387             {
5388                 ind->nsend[zone] = 0;
5389             }
5390             ind->nsend[nzone]     = ind->index.size();
5391             ind->nsend[nzone + 1] = comm->dth[0].nat;
5392             /* Communicate the number of cg's and atoms to receive */
5393             ddSendrecv(dd, dim_ind, dddirBackward,
5394                        ind->nsend, nzone+2,
5395                        ind->nrecv, nzone+2);
5396
5397             if (p > 0)
5398             {
5399                 /* We can receive in place if only the last zone is not empty */
5400                 for (zone = 0; zone < nzone-1; zone++)
5401                 {
5402                     if (ind->nrecv[zone] > 0)
5403                     {
5404                         cd->receiveInPlace = false;
5405                     }
5406                 }
5407             }
5408
5409             int receiveBufferSize = 0;
5410             if (!cd->receiveInPlace)
5411             {
5412                 receiveBufferSize = ind->nrecv[nzone];
5413             }
5414             /* These buffer are actually only needed with in-place */
5415             DDBufferAccess<int>       globalAtomGroupBuffer(comm->intBuffer, receiveBufferSize);
5416             DDBufferAccess<gmx::RVec> rvecBuffer(comm->rvecBuffer, receiveBufferSize);
5417
5418             dd_comm_setup_work_t     &work = comm->dth[0];
5419
5420             /* Make space for the global cg indices */
5421             int numAtomGroupsNew = pos_cg + ind->nrecv[nzone];
5422             dd->globalAtomGroupIndices.resize(numAtomGroupsNew);
5423             /* Communicate the global cg indices */
5424             gmx::ArrayRef<int> integerBufferRef;
5425             if (cd->receiveInPlace)
5426             {
5427                 integerBufferRef = gmx::arrayRefFromArray(dd->globalAtomGroupIndices.data() + pos_cg, ind->nrecv[nzone]);
5428             }
5429             else
5430             {
5431                 integerBufferRef = globalAtomGroupBuffer.buffer;
5432             }
5433             ddSendrecv<int>(dd, dim_ind, dddirBackward,
5434                             work.atomGroupBuffer, integerBufferRef);
5435
5436             /* Make space for cg_cm */
5437             dd_check_alloc_ncg(fr, state, f, pos_cg + ind->nrecv[nzone]);
5438             if (fr->cutoff_scheme == ecutsGROUP)
5439             {
5440                 cg_cm = fr->cg_cm;
5441             }
5442             else
5443             {
5444                 cg_cm = as_rvec_array(state->x.data());
5445             }
5446             /* Communicate cg_cm */
5447             gmx::ArrayRef<gmx::RVec> rvecBufferRef;
5448             if (cd->receiveInPlace)
5449             {
5450                 rvecBufferRef = gmx::arrayRefFromArray(reinterpret_cast<gmx::RVec *>(cg_cm + pos_cg), ind->nrecv[nzone]);
5451             }
5452             else
5453             {
5454                 rvecBufferRef = rvecBuffer.buffer;
5455             }
5456             ddSendrecv<gmx::RVec>(dd, dim_ind, dddirBackward,
5457                                   work.positionBuffer, rvecBufferRef);
5458
5459             /* Make the charge group index */
5460             if (cd->receiveInPlace)
5461             {
5462                 zone = (p == 0 ? 0 : nzone - 1);
5463                 while (zone < nzone)
5464                 {
5465                     for (cg = 0; cg < ind->nrecv[zone]; cg++)
5466                     {
5467                         cg_gl                              = dd->globalAtomGroupIndices[pos_cg];
5468                         fr->cginfo[pos_cg]                 = ddcginfo(cginfo_mb, cg_gl);
5469                         nrcg                               = GET_CGINFO_NATOMS(fr->cginfo[pos_cg]);
5470                         dd->atomGrouping_.appendBlock(nrcg);
5471                         if (bBondComm)
5472                         {
5473                             /* Update the charge group presence,
5474                              * so we can use it in the next pass of the loop.
5475                              */
5476                             comm->bLocalCG[cg_gl] = TRUE;
5477                         }
5478                         pos_cg++;
5479                     }
5480                     if (p == 0)
5481                     {
5482                         comm->zone_ncg1[nzone+zone] = ind->nrecv[zone];
5483                     }
5484                     zone++;
5485                     zone_cg_range[nzone+zone] = pos_cg;
5486                 }
5487             }
5488             else
5489             {
5490                 /* This part of the code is never executed with bBondComm. */
5491                 std::vector<int> &atomGroupsIndex = dd->atomGrouping_.rawIndex();
5492                 atomGroupsIndex.resize(numAtomGroupsNew + 1);
5493
5494                 merge_cg_buffers(nzone, cd, p, zone_cg_range,
5495                                  dd->globalAtomGroupIndices, integerBufferRef.data(),
5496                                  cg_cm, as_rvec_array(rvecBufferRef.data()),
5497                                  atomGroupsIndex,
5498                                  fr->cginfo_mb, fr->cginfo);
5499                 pos_cg += ind->nrecv[nzone];
5500             }
5501             nat_tot += ind->nrecv[nzone+1];
5502         }
5503         if (!cd->receiveInPlace)
5504         {
5505             /* Store the atom block for easy copying of communication buffers */
5506             make_cell2at_index(cd, nzone, zone_cg_range[nzone], dd->atomGrouping());
5507         }
5508         nzone += nzone;
5509     }
5510
5511     comm->atomRanges.setEnd(DDAtomRanges::Type::Zones, nat_tot);
5512
5513     if (!bBondComm)
5514     {
5515         /* We don't need to update cginfo, since that was alrady done above.
5516          * So we pass NULL for the forcerec.
5517          */
5518         dd_set_cginfo(dd->globalAtomGroupIndices,
5519                       dd->ncg_home, dd->globalAtomGroupIndices.size(),
5520                       nullptr, comm->bLocalCG);
5521     }
5522
5523     if (debug)
5524     {
5525         fprintf(debug, "Finished setting up DD communication, zones:");
5526         for (c = 0; c < zones->n; c++)
5527         {
5528             fprintf(debug, " %d", zones->cg_range[c+1]-zones->cg_range[c]);
5529         }
5530         fprintf(debug, "\n");
5531     }
5532 }
5533
5534 static void set_cg_boundaries(gmx_domdec_zones_t *zones)
5535 {
5536     int c;
5537
5538     for (c = 0; c < zones->nizone; c++)
5539     {
5540         zones->izone[c].cg1  = zones->cg_range[c+1];
5541         zones->izone[c].jcg0 = zones->cg_range[zones->izone[c].j0];
5542         zones->izone[c].jcg1 = zones->cg_range[zones->izone[c].j1];
5543     }
5544 }
5545
5546 /* \brief Set zone dimensions for zones \p zone_start to \p zone_end-1
5547  *
5548  * Also sets the atom density for the home zone when \p zone_start=0.
5549  * For this \p numMovedChargeGroupsInHomeZone needs to be passed to tell
5550  * how many charge groups will move but are still part of the current range.
5551  * \todo When converting domdec to use proper classes, all these variables
5552  *       should be private and a method should return the correct count
5553  *       depending on an internal state.
5554  *
5555  * \param[in,out] dd          The domain decomposition struct
5556  * \param[in]     box         The box
5557  * \param[in]     ddbox       The domain decomposition box struct
5558  * \param[in]     zone_start  The start of the zone range to set sizes for
5559  * \param[in]     zone_end    The end of the zone range to set sizes for
5560  * \param[in]     numMovedChargeGroupsInHomeZone  The number of charge groups in the home zone that should moved but are still present in dd->comm->zones.cg_range
5561  */
5562 static void set_zones_size(gmx_domdec_t *dd,
5563                            matrix box, const gmx_ddbox_t *ddbox,
5564                            int zone_start, int zone_end,
5565                            int numMovedChargeGroupsInHomeZone)
5566 {
5567     gmx_domdec_comm_t  *comm;
5568     gmx_domdec_zones_t *zones;
5569     gmx_bool            bDistMB;
5570     int                 z, zi, d, dim;
5571     real                rcs, rcmbs;
5572     int                 i, j;
5573     real                vol;
5574
5575     comm = dd->comm;
5576
5577     zones = &comm->zones;
5578
5579     /* Do we need to determine extra distances for multi-body bondeds? */
5580     bDistMB = (comm->bInterCGMultiBody && isDlbOn(dd->comm) && dd->ndim > 1);
5581
5582     for (z = zone_start; z < zone_end; z++)
5583     {
5584         /* Copy cell limits to zone limits.
5585          * Valid for non-DD dims and non-shifted dims.
5586          */
5587         copy_rvec(comm->cell_x0, zones->size[z].x0);
5588         copy_rvec(comm->cell_x1, zones->size[z].x1);
5589     }
5590
5591     for (d = 0; d < dd->ndim; d++)
5592     {
5593         dim = dd->dim[d];
5594
5595         for (z = 0; z < zones->n; z++)
5596         {
5597             /* With a staggered grid we have different sizes
5598              * for non-shifted dimensions.
5599              */
5600             if (isDlbOn(dd->comm) && zones->shift[z][dim] == 0)
5601             {
5602                 if (d == 1)
5603                 {
5604                     zones->size[z].x0[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min0;
5605                     zones->size[z].x1[dim] = comm->zone_d1[zones->shift[z][dd->dim[d-1]]].max1;
5606                 }
5607                 else if (d == 2)
5608                 {
5609                     zones->size[z].x0[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min0;
5610                     zones->size[z].x1[dim] = comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].max1;
5611                 }
5612             }
5613         }
5614
5615         rcs   = comm->cutoff;
5616         rcmbs = comm->cutoff_mbody;
5617         if (ddbox->tric_dir[dim])
5618         {
5619             rcs   /= ddbox->skew_fac[dim];
5620             rcmbs /= ddbox->skew_fac[dim];
5621         }
5622
5623         /* Set the lower limit for the shifted zone dimensions */
5624         for (z = zone_start; z < zone_end; z++)
5625         {
5626             if (zones->shift[z][dim] > 0)
5627             {
5628                 dim = dd->dim[d];
5629                 if (!isDlbOn(dd->comm) || d == 0)
5630                 {
5631                     zones->size[z].x0[dim] = comm->cell_x1[dim];
5632                     zones->size[z].x1[dim] = comm->cell_x1[dim] + rcs;
5633                 }
5634                 else
5635                 {
5636                     /* Here we take the lower limit of the zone from
5637                      * the lowest domain of the zone below.
5638                      */
5639                     if (z < 4)
5640                     {
5641                         zones->size[z].x0[dim] =
5642                             comm->zone_d1[zones->shift[z][dd->dim[d-1]]].min1;
5643                     }
5644                     else
5645                     {
5646                         if (d == 1)
5647                         {
5648                             zones->size[z].x0[dim] =
5649                                 zones->size[zone_perm[2][z-4]].x0[dim];
5650                         }
5651                         else
5652                         {
5653                             zones->size[z].x0[dim] =
5654                                 comm->zone_d2[zones->shift[z][dd->dim[d-2]]][zones->shift[z][dd->dim[d-1]]].min1;
5655                         }
5656                     }
5657                     /* A temporary limit, is updated below */
5658                     zones->size[z].x1[dim] = zones->size[z].x0[dim];
5659
5660                     if (bDistMB)
5661                     {
5662                         for (zi = 0; zi < zones->nizone; zi++)
5663                         {
5664                             if (zones->shift[zi][dim] == 0)
5665                             {
5666                                 /* This takes the whole zone into account.
5667                                  * With multiple pulses this will lead
5668                                  * to a larger zone then strictly necessary.
5669                                  */
5670                                 zones->size[z].x1[dim] = std::max(zones->size[z].x1[dim],
5671                                                                   zones->size[zi].x1[dim]+rcmbs);
5672                             }
5673                         }
5674                     }
5675                 }
5676             }
5677         }
5678
5679         /* Loop over the i-zones to set the upper limit of each
5680          * j-zone they see.
5681          */
5682         for (zi = 0; zi < zones->nizone; zi++)
5683         {
5684             if (zones->shift[zi][dim] == 0)
5685             {
5686                 /* We should only use zones up to zone_end */
5687                 int jZoneEnd = std::min(zones->izone[zi].j1, zone_end);
5688                 for (z = zones->izone[zi].j0; z < jZoneEnd; z++)
5689                 {
5690                     if (zones->shift[z][dim] > 0)
5691                     {
5692                         zones->size[z].x1[dim] = std::max(zones->size[z].x1[dim],
5693                                                           zones->size[zi].x1[dim]+rcs);
5694                     }
5695                 }
5696             }
5697         }
5698     }
5699
5700     for (z = zone_start; z < zone_end; z++)
5701     {
5702         /* Initialization only required to keep the compiler happy */
5703         rvec corner_min = {0, 0, 0}, corner_max = {0, 0, 0}, corner;
5704         int  nc, c;
5705
5706         /* To determine the bounding box for a zone we need to find
5707          * the extreme corners of 4, 2 or 1 corners.
5708          */
5709         nc = 1 << (ddbox->nboundeddim - 1);
5710
5711         for (c = 0; c < nc; c++)
5712         {
5713             /* Set up a zone corner at x=0, ignoring trilinic couplings */
5714             corner[XX] = 0;
5715             if ((c & 1) == 0)
5716             {
5717                 corner[YY] = zones->size[z].x0[YY];
5718             }
5719             else
5720             {
5721                 corner[YY] = zones->size[z].x1[YY];
5722             }
5723             if ((c & 2) == 0)
5724             {
5725                 corner[ZZ] = zones->size[z].x0[ZZ];
5726             }
5727             else
5728             {
5729                 corner[ZZ] = zones->size[z].x1[ZZ];
5730             }
5731             if (dd->ndim == 1 && dd->dim[0] < ZZ && ZZ < dd->npbcdim &&
5732                 box[ZZ][1 - dd->dim[0]] != 0)
5733             {
5734                 /* With 1D domain decomposition the cg's are not in
5735                  * the triclinic box, but triclinic x-y and rectangular y/x-z.
5736                  * Shift the corner of the z-vector back to along the box
5737                  * vector of dimension d, so it will later end up at 0 along d.
5738                  * This can affect the location of this corner along dd->dim[0]
5739                  * through the matrix operation below if box[d][dd->dim[0]]!=0.
5740                  */
5741                 int d = 1 - dd->dim[0];
5742
5743                 corner[d] -= corner[ZZ]*box[ZZ][d]/box[ZZ][ZZ];
5744             }
5745             /* Apply the triclinic couplings */
5746             assert(ddbox->npbcdim <= DIM);
5747             for (i = YY; i < ddbox->npbcdim; i++)
5748             {
5749                 for (j = XX; j < i; j++)
5750                 {
5751                     corner[j] += corner[i]*box[i][j]/box[i][i];
5752                 }
5753             }
5754             if (c == 0)
5755             {
5756                 copy_rvec(corner, corner_min);
5757                 copy_rvec(corner, corner_max);
5758             }
5759             else
5760             {
5761                 for (i = 0; i < DIM; i++)
5762                 {
5763                     corner_min[i] = std::min(corner_min[i], corner[i]);
5764                     corner_max[i] = std::max(corner_max[i], corner[i]);
5765                 }
5766             }
5767         }
5768         /* Copy the extreme cornes without offset along x */
5769         for (i = 0; i < DIM; i++)
5770         {
5771             zones->size[z].bb_x0[i] = corner_min[i];
5772             zones->size[z].bb_x1[i] = corner_max[i];
5773         }
5774         /* Add the offset along x */
5775         zones->size[z].bb_x0[XX] += zones->size[z].x0[XX];
5776         zones->size[z].bb_x1[XX] += zones->size[z].x1[XX];
5777     }
5778
5779     if (zone_start == 0)
5780     {
5781         vol = 1;
5782         for (dim = 0; dim < DIM; dim++)
5783         {
5784             vol *= zones->size[0].x1[dim] - zones->size[0].x0[dim];
5785         }
5786         zones->dens_zone0 = (zones->cg_range[1] - zones->cg_range[0] - numMovedChargeGroupsInHomeZone)/vol;
5787     }
5788
5789     if (debug)
5790     {
5791         for (z = zone_start; z < zone_end; z++)
5792         {
5793             fprintf(debug, "zone %d    %6.3f - %6.3f  %6.3f - %6.3f  %6.3f - %6.3f\n",
5794                     z,
5795                     zones->size[z].x0[XX], zones->size[z].x1[XX],
5796                     zones->size[z].x0[YY], zones->size[z].x1[YY],
5797                     zones->size[z].x0[ZZ], zones->size[z].x1[ZZ]);
5798             fprintf(debug, "zone %d bb %6.3f - %6.3f  %6.3f - %6.3f  %6.3f - %6.3f\n",
5799                     z,
5800                     zones->size[z].bb_x0[XX], zones->size[z].bb_x1[XX],
5801                     zones->size[z].bb_x0[YY], zones->size[z].bb_x1[YY],
5802                     zones->size[z].bb_x0[ZZ], zones->size[z].bb_x1[ZZ]);
5803         }
5804     }
5805 }
5806
5807 static int comp_cgsort(const void *a, const void *b)
5808 {
5809     int           comp;
5810
5811     gmx_cgsort_t *cga, *cgb;
5812     cga = (gmx_cgsort_t *)a;
5813     cgb = (gmx_cgsort_t *)b;
5814
5815     comp = cga->nsc - cgb->nsc;
5816     if (comp == 0)
5817     {
5818         comp = cga->ind_gl - cgb->ind_gl;
5819     }
5820
5821     return comp;
5822 }
5823
5824 /* Order data in \p dataToSort according to \p sort
5825  *
5826  * Note: both buffers should have at least \p sort.size() elements.
5827  */
5828 template <typename T>
5829 static void
5830 orderVector(gmx::ArrayRef<const gmx_cgsort_t>  sort,
5831             gmx::ArrayRef<T>                   dataToSort,
5832             gmx::ArrayRef<T>                   sortBuffer)
5833 {
5834     GMX_ASSERT(dataToSort.size() >= sort.size(), "The vector needs to be sufficiently large");
5835     GMX_ASSERT(sortBuffer.size() >= sort.size(), "The sorting buffer needs to be sufficiently large");
5836
5837     /* Order the data into the temporary buffer */
5838     size_t i = 0;
5839     for (const gmx_cgsort_t &entry : sort)
5840     {
5841         sortBuffer[i++] = dataToSort[entry.ind];
5842     }
5843
5844     /* Copy back to the original array */
5845     std::copy(sortBuffer.begin(), sortBuffer.begin() + sort.size(),
5846               dataToSort.begin());
5847 }
5848
5849 /* Order data in \p dataToSort according to \p sort
5850  *
5851  * Note: \p vectorToSort should have at least \p sort.size() elements,
5852  *       \p workVector is resized when it is too small.
5853  */
5854 template <typename T>
5855 static void
5856 orderVector(gmx::ArrayRef<const gmx_cgsort_t>  sort,
5857             gmx::ArrayRef<T>                   vectorToSort,
5858             std::vector<T>                    *workVector)
5859 {
5860     if (gmx::index(workVector->size()) < sort.size())
5861     {
5862         workVector->resize(sort.size());
5863     }
5864     orderVector<T>(sort, vectorToSort, *workVector);
5865 }
5866
5867 static void order_vec_atom(const gmx::RangePartitioning      *atomGroups,
5868                            gmx::ArrayRef<const gmx_cgsort_t>  sort,
5869                            gmx::ArrayRef<gmx::RVec>           v,
5870                            gmx::ArrayRef<gmx::RVec>           buf)
5871 {
5872     if (atomGroups == nullptr)
5873     {
5874         /* Avoid the useless loop of the atoms within a cg */
5875         orderVector(sort, v, buf);
5876
5877         return;
5878     }
5879
5880     /* Order the data */
5881     int a = 0;
5882     for (const gmx_cgsort_t &entry : sort)
5883     {
5884         for (int i : atomGroups->block(entry.ind))
5885         {
5886             copy_rvec(v[i], buf[a]);
5887             a++;
5888         }
5889     }
5890     int atot = a;
5891
5892     /* Copy back to the original array */
5893     for (int a = 0; a < atot; a++)
5894     {
5895         copy_rvec(buf[a], v[a]);
5896     }
5897 }
5898
5899 /* Returns whether a < b */
5900 static bool compareCgsort(const gmx_cgsort_t &a,
5901                           const gmx_cgsort_t &b)
5902 {
5903     return (a.nsc < b.nsc ||
5904             (a.nsc == b.nsc && a.ind_gl < b.ind_gl));
5905 }
5906
5907 static void orderedSort(gmx::ArrayRef<const gmx_cgsort_t>  stationary,
5908                         gmx::ArrayRef<gmx_cgsort_t>        moved,
5909                         std::vector<gmx_cgsort_t>         *sort1)
5910 {
5911     /* The new indices are not very ordered, so we qsort them */
5912     gmx_qsort_threadsafe(moved.data(), moved.size(), sizeof(moved[0]), comp_cgsort);
5913
5914     /* stationary is already ordered, so now we can merge the two arrays */
5915     sort1->resize(stationary.size() + moved.size());
5916     std::merge(stationary.begin(), stationary.end(),
5917                moved.begin(), moved.end(),
5918                sort1->begin(),
5919                compareCgsort);
5920 }
5921
5922 /* Set the sorting order for systems with charge groups, returned in sort->sort.
5923  * The order is according to the global charge group index.
5924  * This adds and removes charge groups that moved between domains.
5925  */
5926 static void dd_sort_order(const gmx_domdec_t *dd,
5927                           const t_forcerec   *fr,
5928                           int                 ncg_home_old,
5929                           gmx_domdec_sort_t  *sort)
5930 {
5931     const int *a          = fr->ns->grid->cell_index;
5932
5933     const int  movedValue = NSGRID_SIGNAL_MOVED_FAC*fr->ns->grid->ncells;
5934
5935     if (ncg_home_old >= 0)
5936     {
5937         std::vector<gmx_cgsort_t> &stationary = sort->stationary;
5938         std::vector<gmx_cgsort_t> &moved      = sort->moved;
5939
5940         /* The charge groups that remained in the same ns grid cell
5941          * are completely ordered. So we can sort efficiently by sorting
5942          * the charge groups that did move into the stationary list.
5943          * Note: push_back() seems to be slightly slower than direct access.
5944          */
5945         stationary.clear();
5946         moved.clear();
5947         for (int i = 0; i < dd->ncg_home; i++)
5948         {
5949             /* Check if this cg did not move to another node */
5950             if (a[i] < movedValue)
5951             {
5952                 gmx_cgsort_t entry;
5953                 entry.nsc    = a[i];
5954                 entry.ind_gl = dd->globalAtomGroupIndices[i];
5955                 entry.ind    = i;
5956
5957                 if (i >= ncg_home_old || a[i] != sort->sorted[i].nsc)
5958                 {
5959                     /* This cg is new on this node or moved ns grid cell */
5960                     moved.push_back(entry);
5961                 }
5962                 else
5963                 {
5964                     /* This cg did not move */
5965                     stationary.push_back(entry);
5966                 }
5967             }
5968         }
5969
5970         if (debug)
5971         {
5972             fprintf(debug, "ordered sort cgs: stationary %zu moved %zu\n",
5973                     stationary.size(), moved.size());
5974         }
5975         /* Sort efficiently */
5976         orderedSort(stationary, moved, &sort->sorted);
5977     }
5978     else
5979     {
5980         std::vector<gmx_cgsort_t> &cgsort   = sort->sorted;
5981         cgsort.clear();
5982         cgsort.reserve(dd->ncg_home);
5983         int                        numCGNew = 0;
5984         for (int i = 0; i < dd->ncg_home; i++)
5985         {
5986             /* Sort on the ns grid cell indices
5987              * and the global topology index
5988              */
5989             gmx_cgsort_t entry;
5990             entry.nsc    = a[i];
5991             entry.ind_gl = dd->globalAtomGroupIndices[i];
5992             entry.ind    = i;
5993             cgsort.push_back(entry);
5994             if (cgsort[i].nsc < movedValue)
5995             {
5996                 numCGNew++;
5997             }
5998         }
5999         if (debug)
6000         {
6001             fprintf(debug, "qsort cgs: %d new home %d\n", dd->ncg_home, numCGNew);
6002         }
6003         /* Determine the order of the charge groups using qsort */
6004         gmx_qsort_threadsafe(cgsort.data(), dd->ncg_home, sizeof(cgsort[0]), comp_cgsort);
6005
6006         /* Remove the charge groups which are no longer at home here */
6007         cgsort.resize(numCGNew);
6008     }
6009 }
6010
6011 /* Returns the sorting order for atoms based on the nbnxn grid order in sort */
6012 static void dd_sort_order_nbnxn(const t_forcerec          *fr,
6013                                 std::vector<gmx_cgsort_t> *sort)
6014 {
6015     gmx::ArrayRef<const int> atomOrder = nbnxn_get_atomorder(fr->nbv->nbs.get());
6016
6017     /* Using push_back() instead of this resize results in much slower code */
6018     sort->resize(atomOrder.size());
6019     gmx::ArrayRef<gmx_cgsort_t> buffer    = *sort;
6020     size_t                      numSorted = 0;
6021     for (int i : atomOrder)
6022     {
6023         if (i >= 0)
6024         {
6025             /* The values of nsc and ind_gl are not used in this case */
6026             buffer[numSorted++].ind = i;
6027         }
6028     }
6029     sort->resize(numSorted);
6030 }
6031
6032 static void dd_sort_state(gmx_domdec_t *dd, rvec *cgcm, t_forcerec *fr, t_state *state,
6033                           int ncg_home_old)
6034 {
6035     gmx_domdec_sort_t *sort = dd->comm->sort.get();
6036
6037     switch (fr->cutoff_scheme)
6038     {
6039         case ecutsGROUP:
6040             dd_sort_order(dd, fr, ncg_home_old, sort);
6041             break;
6042         case ecutsVERLET:
6043             dd_sort_order_nbnxn(fr, &sort->sorted);
6044             break;
6045         default:
6046             gmx_incons("unimplemented");
6047     }
6048
6049     const gmx::RangePartitioning &atomGrouping = dd->atomGrouping();
6050
6051     /* We alloc with the old size, since cgindex is still old */
6052     GMX_ASSERT(atomGrouping.numBlocks() == dd->ncg_home, "atomGroups and dd should be consistent");
6053     DDBufferAccess<gmx::RVec>     rvecBuffer(dd->comm->rvecBuffer, atomGrouping.fullRange().end());
6054
6055     const gmx::RangePartitioning *atomGroupsPtr = (dd->comm->bCGs ? &atomGrouping : nullptr);
6056
6057     /* Set the new home atom/charge group count */
6058     dd->ncg_home = sort->sorted.size();
6059     if (debug)
6060     {
6061         fprintf(debug, "Set the new home charge group count to %d\n",
6062                 dd->ncg_home);
6063     }
6064
6065     /* Reorder the state */
6066     gmx::ArrayRef<const gmx_cgsort_t> cgsort = sort->sorted;
6067     GMX_RELEASE_ASSERT(cgsort.size() == dd->ncg_home, "We should sort all the home atom groups");
6068
6069     if (state->flags & (1 << estX))
6070     {
6071         order_vec_atom(atomGroupsPtr, cgsort, state->x, rvecBuffer.buffer);
6072     }
6073     if (state->flags & (1 << estV))
6074     {
6075         order_vec_atom(atomGroupsPtr, cgsort, state->v, rvecBuffer.buffer);
6076     }
6077     if (state->flags & (1 << estCGP))
6078     {
6079         order_vec_atom(atomGroupsPtr, cgsort, state->cg_p, rvecBuffer.buffer);
6080     }
6081
6082     if (fr->cutoff_scheme == ecutsGROUP)
6083     {
6084         /* Reorder cgcm */
6085         gmx::ArrayRef<gmx::RVec> cgcmRef = gmx::arrayRefFromArray(reinterpret_cast<gmx::RVec *>(cgcm[0]), cgsort.size());
6086         orderVector(cgsort, cgcmRef, rvecBuffer.buffer);
6087     }
6088
6089     /* Reorder the global cg index */
6090     orderVector<int>(cgsort, dd->globalAtomGroupIndices, &sort->intBuffer);
6091     /* Reorder the cginfo */
6092     orderVector<int>(cgsort, gmx::arrayRefFromArray(fr->cginfo, cgsort.size()), &sort->intBuffer);
6093     /* Rebuild the local cg index */
6094     if (dd->comm->bCGs)
6095     {
6096         /* We make a new, ordered atomGroups object and assign it to
6097          * the old one. This causes some allocation overhead, but saves
6098          * a copy back of the whole index.
6099          */
6100         gmx::RangePartitioning ordered;
6101         for (const gmx_cgsort_t &entry : cgsort)
6102         {
6103             ordered.appendBlock(atomGrouping.block(entry.ind).size());
6104         }
6105         dd->atomGrouping_ = ordered;
6106     }
6107     else
6108     {
6109         dd->atomGrouping_.setAllBlocksSizeOne(dd->ncg_home);
6110     }
6111     /* Set the home atom number */
6112     dd->comm->atomRanges.setEnd(DDAtomRanges::Type::Home, dd->atomGrouping().fullRange().end());
6113
6114     if (fr->cutoff_scheme == ecutsVERLET)
6115     {
6116         /* The atoms are now exactly in grid order, update the grid order */
6117         nbnxn_set_atomorder(fr->nbv->nbs.get());
6118     }
6119     else
6120     {
6121         /* Copy the sorted ns cell indices back to the ns grid struct */
6122         for (gmx::index i = 0; i < cgsort.size(); i++)
6123         {
6124             fr->ns->grid->cell_index[i] = cgsort[i].nsc;
6125         }
6126         fr->ns->grid->nr = cgsort.size();
6127     }
6128 }
6129
6130 static void add_dd_statistics(gmx_domdec_t *dd)
6131 {
6132     gmx_domdec_comm_t *comm = dd->comm;
6133
6134     for (int i = 0; i < static_cast<int>(DDAtomRanges::Type::Number); i++)
6135     {
6136         auto range = static_cast<DDAtomRanges::Type>(i);
6137         comm->sum_nat[i] +=
6138             comm->atomRanges.end(range) - comm->atomRanges.start(range);
6139     }
6140     comm->ndecomp++;
6141 }
6142
6143 void reset_dd_statistics_counters(gmx_domdec_t *dd)
6144 {
6145     gmx_domdec_comm_t *comm = dd->comm;
6146
6147     /* Reset all the statistics and counters for total run counting */
6148     for (int i = 0; i < static_cast<int>(DDAtomRanges::Type::Number); i++)
6149     {
6150         comm->sum_nat[i] = 0;
6151     }
6152     comm->ndecomp   = 0;
6153     comm->nload     = 0;
6154     comm->load_step = 0;
6155     comm->load_sum  = 0;
6156     comm->load_max  = 0;
6157     clear_ivec(comm->load_lim);
6158     comm->load_mdf = 0;
6159     comm->load_pme = 0;
6160 }
6161
6162 void print_dd_statistics(const t_commrec *cr, const t_inputrec *ir, FILE *fplog)
6163 {
6164     gmx_domdec_comm_t *comm      = cr->dd->comm;
6165
6166     const int          numRanges = static_cast<int>(DDAtomRanges::Type::Number);
6167     gmx_sumd(numRanges, comm->sum_nat, cr);
6168
6169     if (fplog == nullptr)
6170     {
6171         return;
6172     }
6173
6174     fprintf(fplog, "\n    D O M A I N   D E C O M P O S I T I O N   S T A T I S T I C S\n\n");
6175
6176     for (int i = static_cast<int>(DDAtomRanges::Type::Zones); i < numRanges; i++)
6177     {
6178         auto   range = static_cast<DDAtomRanges::Type>(i);
6179         double av    = comm->sum_nat[i]/comm->ndecomp;
6180         switch (range)
6181         {
6182             case DDAtomRanges::Type::Zones:
6183                 fprintf(fplog,
6184                         " av. #atoms communicated per step for force:  %d x %.1f\n",
6185                         2, av);
6186                 break;
6187             case DDAtomRanges::Type::Vsites:
6188                 if (cr->dd->vsite_comm)
6189                 {
6190                     fprintf(fplog,
6191                             " av. #atoms communicated per step for vsites: %d x %.1f\n",
6192                             (EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD) ? 3 : 2,
6193                             av);
6194                 }
6195                 break;
6196             case DDAtomRanges::Type::Constraints:
6197                 if (cr->dd->constraint_comm)
6198                 {
6199                     fprintf(fplog,
6200                             " av. #atoms communicated per step for LINCS:  %d x %.1f\n",
6201                             1 + ir->nLincsIter, av);
6202                 }
6203                 break;
6204             default:
6205                 gmx_incons(" Unknown type for DD statistics");
6206         }
6207     }
6208     fprintf(fplog, "\n");
6209
6210     if (comm->bRecordLoad && EI_DYNAMICS(ir->eI))
6211     {
6212         print_dd_load_av(fplog, cr->dd);
6213     }
6214 }
6215
6216 void dd_partition_system(FILE                *fplog,
6217                          int64_t              step,
6218                          const t_commrec     *cr,
6219                          gmx_bool             bMasterState,
6220                          int                  nstglobalcomm,
6221                          t_state             *state_global,
6222                          const gmx_mtop_t    *top_global,
6223                          const t_inputrec    *ir,
6224                          gmx_enfrot          *enforcedRotation,
6225                          t_state             *state_local,
6226                          PaddedRVecVector    *f,
6227                          gmx::MDAtoms        *mdAtoms,
6228                          gmx_localtop_t      *top_local,
6229                          t_forcerec          *fr,
6230                          gmx_vsite_t         *vsite,
6231                          gmx::Constraints    *constr,
6232                          t_nrnb              *nrnb,
6233                          gmx_wallcycle       *wcycle,
6234                          gmx_bool             bVerbose)
6235 {
6236     gmx_domdec_t      *dd;
6237     gmx_domdec_comm_t *comm;
6238     gmx_ddbox_t        ddbox = {0};
6239     t_block           *cgs_gl;
6240     int64_t            step_pcoupl;
6241     rvec               cell_ns_x0, cell_ns_x1;
6242     int                ncgindex_set, ncg_home_old = -1, ncg_moved, nat_f_novirsum;
6243     gmx_bool           bBoxChanged, bNStGlobalComm, bDoDLB, bCheckWhetherToTurnDlbOn, bLogLoad;
6244     gmx_bool           bRedist, bSortCG, bResortAll;
6245     ivec               ncells_old = {0, 0, 0}, ncells_new = {0, 0, 0}, np;
6246     real               grid_density;
6247     char               sbuf[22];
6248
6249     wallcycle_start(wcycle, ewcDOMDEC);
6250
6251     dd   = cr->dd;
6252     comm = dd->comm;
6253
6254     // TODO if the update code becomes accessible here, use
6255     // upd->deform for this logic.
6256     bBoxChanged = (bMasterState || inputrecDeform(ir));
6257     if (ir->epc != epcNO)
6258     {
6259         /* With nstpcouple > 1 pressure coupling happens.
6260          * one step after calculating the pressure.
6261          * Box scaling happens at the end of the MD step,
6262          * after the DD partitioning.
6263          * We therefore have to do DLB in the first partitioning
6264          * after an MD step where P-coupling occurred.
6265          * We need to determine the last step in which p-coupling occurred.
6266          * MRS -- need to validate this for vv?
6267          */
6268         int n = ir->nstpcouple;
6269         if (n == 1)
6270         {
6271             step_pcoupl = step - 1;
6272         }
6273         else
6274         {
6275             step_pcoupl = ((step - 1)/n)*n + 1;
6276         }
6277         if (step_pcoupl >= comm->partition_step)
6278         {
6279             bBoxChanged = TRUE;
6280         }
6281     }
6282
6283     bNStGlobalComm = (step % nstglobalcomm == 0);
6284
6285     if (!isDlbOn(comm))
6286     {
6287         bDoDLB = FALSE;
6288     }
6289     else
6290     {
6291         /* Should we do dynamic load balacing this step?
6292          * Since it requires (possibly expensive) global communication,
6293          * we might want to do DLB less frequently.
6294          */
6295         if (bBoxChanged || ir->epc != epcNO)
6296         {
6297             bDoDLB = bBoxChanged;
6298         }
6299         else
6300         {
6301             bDoDLB = bNStGlobalComm;
6302         }
6303     }
6304
6305     /* Check if we have recorded loads on the nodes */
6306     if (comm->bRecordLoad && dd_load_count(comm) > 0)
6307     {
6308         bCheckWhetherToTurnDlbOn = dd_dlb_get_should_check_whether_to_turn_dlb_on(dd);
6309
6310         /* Print load every nstlog, first and last step to the log file */
6311         bLogLoad = ((ir->nstlog > 0 && step % ir->nstlog == 0) ||
6312                     comm->n_load_collect == 0 ||
6313                     (ir->nsteps >= 0 &&
6314                      (step + ir->nstlist > ir->init_step + ir->nsteps)));
6315
6316         /* Avoid extra communication due to verbose screen output
6317          * when nstglobalcomm is set.
6318          */
6319         if (bDoDLB || bLogLoad || bCheckWhetherToTurnDlbOn ||
6320             (bVerbose && (ir->nstlist == 0 || nstglobalcomm <= ir->nstlist)))
6321         {
6322             get_load_distribution(dd, wcycle);
6323             if (DDMASTER(dd))
6324             {
6325                 if (bLogLoad)
6326                 {
6327                     dd_print_load(fplog, dd, step-1);
6328                 }
6329                 if (bVerbose)
6330                 {
6331                     dd_print_load_verbose(dd);
6332                 }
6333             }
6334             comm->n_load_collect++;
6335
6336             if (isDlbOn(comm))
6337             {
6338                 if (DDMASTER(dd))
6339                 {
6340                     /* Add the measured cycles to the running average */
6341                     const float averageFactor        = 0.1f;
6342                     comm->cyclesPerStepDlbExpAverage =
6343                         (1 - averageFactor)*comm->cyclesPerStepDlbExpAverage +
6344                         averageFactor*comm->cycl[ddCyclStep]/comm->cycl_n[ddCyclStep];
6345                 }
6346                 if (comm->dlbState == edlbsOnCanTurnOff &&
6347                     dd->comm->n_load_have % c_checkTurnDlbOffInterval == c_checkTurnDlbOffInterval - 1)
6348                 {
6349                     gmx_bool turnOffDlb;
6350                     if (DDMASTER(dd))
6351                     {
6352                         /* If the running averaged cycles with DLB are more
6353                          * than before we turned on DLB, turn off DLB.
6354                          * We will again run and check the cycles without DLB
6355                          * and we can then decide if to turn off DLB forever.
6356                          */
6357                         turnOffDlb = (comm->cyclesPerStepDlbExpAverage >
6358                                       comm->cyclesPerStepBeforeDLB);
6359                     }
6360                     dd_bcast(dd, sizeof(turnOffDlb), &turnOffDlb);
6361                     if (turnOffDlb)
6362                     {
6363                         /* To turn off DLB, we need to redistribute the atoms */
6364                         dd_collect_state(dd, state_local, state_global);
6365                         bMasterState = TRUE;
6366                         turn_off_dlb(fplog, cr, step);
6367                     }
6368                 }
6369             }
6370             else if (bCheckWhetherToTurnDlbOn)
6371             {
6372                 gmx_bool turnOffDlbForever = FALSE;
6373                 gmx_bool turnOnDlb         = FALSE;
6374
6375                 /* Since the timings are node dependent, the master decides */
6376                 if (DDMASTER(dd))
6377                 {
6378                     /* If we recently turned off DLB, we want to check if
6379                      * performance is better without DLB. We want to do this
6380                      * ASAP to minimize the chance that external factors
6381                      * slowed down the DLB step are gone here and we
6382                      * incorrectly conclude that DLB was causing the slowdown.
6383                      * So we measure one nstlist block, no running average.
6384                      */
6385                     if (comm->haveTurnedOffDlb &&
6386                         comm->cycl[ddCyclStep]/comm->cycl_n[ddCyclStep] <
6387                         comm->cyclesPerStepDlbExpAverage)
6388                     {
6389                         /* After turning off DLB we ran nstlist steps in fewer
6390                          * cycles than with DLB. This likely means that DLB
6391                          * in not benefical, but this could be due to a one
6392                          * time unlucky fluctuation, so we require two such
6393                          * observations in close succession to turn off DLB
6394                          * forever.
6395                          */
6396                         if (comm->dlbSlowerPartitioningCount > 0 &&
6397                             dd->ddp_count < comm->dlbSlowerPartitioningCount + 10*c_checkTurnDlbOnInterval)
6398                         {
6399                             turnOffDlbForever = TRUE;
6400                         }
6401                         comm->haveTurnedOffDlb           = false;
6402                         /* Register when we last measured DLB slowdown */
6403                         comm->dlbSlowerPartitioningCount = dd->ddp_count;
6404                     }
6405                     else
6406                     {
6407                         /* Here we check if the max PME rank load is more than 0.98
6408                          * the max PP force load. If so, PP DLB will not help,
6409                          * since we are (almost) limited by PME. Furthermore,
6410                          * DLB will cause a significant extra x/f redistribution
6411                          * cost on the PME ranks, which will then surely result
6412                          * in lower total performance.
6413                          */
6414                         if (cr->npmenodes > 0 &&
6415                             dd_pme_f_ratio(dd) > 1 - DD_PERF_LOSS_DLB_ON)
6416                         {
6417                             turnOnDlb = FALSE;
6418                         }
6419                         else
6420                         {
6421                             turnOnDlb = (dd_force_imb_perf_loss(dd) >= DD_PERF_LOSS_DLB_ON);
6422                         }
6423                     }
6424                 }
6425                 struct
6426                 {
6427                     gmx_bool turnOffDlbForever;
6428                     gmx_bool turnOnDlb;
6429                 }
6430                 bools {
6431                     turnOffDlbForever, turnOnDlb
6432                 };
6433                 dd_bcast(dd, sizeof(bools), &bools);
6434                 if (bools.turnOffDlbForever)
6435                 {
6436                     turn_off_dlb_forever(fplog, cr, step);
6437                 }
6438                 else if (bools.turnOnDlb)
6439                 {
6440                     turn_on_dlb(fplog, cr, step);
6441                     bDoDLB = TRUE;
6442                 }
6443             }
6444         }
6445         comm->n_load_have++;
6446     }
6447
6448     cgs_gl = &comm->cgs_gl;
6449
6450     bRedist = FALSE;
6451     if (bMasterState)
6452     {
6453         /* Clear the old state */
6454         clearDDStateIndices(dd, 0, 0);
6455         ncgindex_set = 0;
6456
6457         auto xGlobal = positionsFromStatePointer(state_global);
6458
6459         set_ddbox(dd, true, ir,
6460                   DDMASTER(dd) ? state_global->box : nullptr,
6461                   true, xGlobal,
6462                   &ddbox);
6463
6464         distributeState(fplog, dd, state_global, ddbox, state_local, f);
6465
6466         dd_make_local_cgs(dd, &top_local->cgs);
6467
6468         /* Ensure that we have space for the new distribution */
6469         dd_check_alloc_ncg(fr, state_local, f, dd->ncg_home);
6470
6471         if (fr->cutoff_scheme == ecutsGROUP)
6472         {
6473             calc_cgcm(fplog, 0, dd->ncg_home,
6474                       &top_local->cgs, as_rvec_array(state_local->x.data()), fr->cg_cm);
6475         }
6476
6477         inc_nrnb(nrnb, eNR_CGCM, comm->atomRanges.numHomeAtoms());
6478
6479         dd_set_cginfo(dd->globalAtomGroupIndices, 0, dd->ncg_home, fr, comm->bLocalCG);
6480     }
6481     else if (state_local->ddp_count != dd->ddp_count)
6482     {
6483         if (state_local->ddp_count > dd->ddp_count)
6484         {
6485             gmx_fatal(FARGS, "Internal inconsistency state_local->ddp_count (%d) > dd->ddp_count (%ld)", state_local->ddp_count, dd->ddp_count);
6486         }
6487
6488         if (state_local->ddp_count_cg_gl != state_local->ddp_count)
6489         {
6490             gmx_fatal(FARGS, "Internal inconsistency state_local->ddp_count_cg_gl (%d) != state_local->ddp_count (%d)", state_local->ddp_count_cg_gl, state_local->ddp_count);
6491         }
6492
6493         /* Clear the old state */
6494         clearDDStateIndices(dd, 0, 0);
6495
6496         /* Restore the atom group indices from state_local */
6497         restoreAtomGroups(dd, cgs_gl->index, state_local);
6498         make_dd_indices(dd, cgs_gl->index, 0);
6499         ncgindex_set = dd->ncg_home;
6500
6501         if (fr->cutoff_scheme == ecutsGROUP)
6502         {
6503             /* Redetermine the cg COMs */
6504             calc_cgcm(fplog, 0, dd->ncg_home,
6505                       &top_local->cgs, as_rvec_array(state_local->x.data()), fr->cg_cm);
6506         }
6507
6508         inc_nrnb(nrnb, eNR_CGCM, comm->atomRanges.numHomeAtoms());
6509
6510         dd_set_cginfo(dd->globalAtomGroupIndices, 0, dd->ncg_home, fr, comm->bLocalCG);
6511
6512         set_ddbox(dd, bMasterState, ir, state_local->box,
6513                   true, state_local->x, &ddbox);
6514
6515         bRedist = isDlbOn(comm);
6516     }
6517     else
6518     {
6519         /* We have the full state, only redistribute the cgs */
6520
6521         /* Clear the non-home indices */
6522         clearDDStateIndices(dd, dd->ncg_home, comm->atomRanges.numHomeAtoms());
6523         ncgindex_set = 0;
6524
6525         /* Avoid global communication for dim's without pbc and -gcom */
6526         if (!bNStGlobalComm)
6527         {
6528             copy_rvec(comm->box0, ddbox.box0    );
6529             copy_rvec(comm->box_size, ddbox.box_size);
6530         }
6531         set_ddbox(dd, bMasterState, ir, state_local->box,
6532                   bNStGlobalComm, state_local->x, &ddbox);
6533
6534         bBoxChanged = TRUE;
6535         bRedist     = TRUE;
6536     }
6537     /* For dim's without pbc and -gcom */
6538     copy_rvec(ddbox.box0, comm->box0    );
6539     copy_rvec(ddbox.box_size, comm->box_size);
6540
6541     set_dd_cell_sizes(dd, &ddbox, dynamic_dd_box(&ddbox, ir), bMasterState, bDoDLB,
6542                       step, wcycle);
6543
6544     if (comm->nstDDDumpGrid > 0 && step % comm->nstDDDumpGrid == 0)
6545     {
6546         write_dd_grid_pdb("dd_grid", step, dd, state_local->box, &ddbox);
6547     }
6548
6549     /* Check if we should sort the charge groups */
6550     bSortCG = (bMasterState || bRedist);
6551
6552     ncg_home_old = dd->ncg_home;
6553
6554     /* When repartitioning we mark charge groups that will move to neighboring
6555      * DD cells, but we do not move them right away for performance reasons.
6556      * Thus we need to keep track of how many charge groups will move for
6557      * obtaining correct local charge group / atom counts.
6558      */
6559     ncg_moved = 0;
6560     if (bRedist)
6561     {
6562         wallcycle_sub_start(wcycle, ewcsDD_REDIST);
6563
6564         dd_redistribute_cg(fplog, step, dd, ddbox.tric_dir,
6565                            state_local, f, fr,
6566                            !bSortCG, nrnb, &ncgindex_set, &ncg_moved);
6567
6568         wallcycle_sub_stop(wcycle, ewcsDD_REDIST);
6569     }
6570
6571     get_nsgrid_boundaries(ddbox.nboundeddim, state_local->box,
6572                           dd, &ddbox,
6573                           &comm->cell_x0, &comm->cell_x1,
6574                           dd->ncg_home, fr->cg_cm,
6575                           cell_ns_x0, cell_ns_x1, &grid_density);
6576
6577     if (bBoxChanged)
6578     {
6579         comm_dd_ns_cell_sizes(dd, &ddbox, cell_ns_x0, cell_ns_x1, step);
6580     }
6581
6582     switch (fr->cutoff_scheme)
6583     {
6584         case ecutsGROUP:
6585             copy_ivec(fr->ns->grid->n, ncells_old);
6586             grid_first(fplog, fr->ns->grid, dd, &ddbox,
6587                        state_local->box, cell_ns_x0, cell_ns_x1,
6588                        fr->rlist, grid_density);
6589             break;
6590         case ecutsVERLET:
6591             nbnxn_get_ncells(fr->nbv->nbs.get(), &ncells_old[XX], &ncells_old[YY]);
6592             break;
6593         default:
6594             gmx_incons("unimplemented");
6595     }
6596     /* We need to store tric_dir for dd_get_ns_ranges called from ns.c */
6597     copy_ivec(ddbox.tric_dir, comm->tric_dir);
6598
6599     if (bSortCG)
6600     {
6601         wallcycle_sub_start(wcycle, ewcsDD_GRID);
6602
6603         /* Sort the state on charge group position.
6604          * This enables exact restarts from this step.
6605          * It also improves performance by about 15% with larger numbers
6606          * of atoms per node.
6607          */
6608
6609         /* Fill the ns grid with the home cell,
6610          * so we can sort with the indices.
6611          */
6612         set_zones_ncg_home(dd);
6613
6614         switch (fr->cutoff_scheme)
6615         {
6616             case ecutsVERLET:
6617                 set_zones_size(dd, state_local->box, &ddbox, 0, 1, ncg_moved);
6618
6619                 nbnxn_put_on_grid(fr->nbv->nbs.get(), fr->ePBC, state_local->box,
6620                                   0,
6621                                   comm->zones.size[0].bb_x0,
6622                                   comm->zones.size[0].bb_x1,
6623                                   0, dd->ncg_home,
6624                                   comm->zones.dens_zone0,
6625                                   fr->cginfo,
6626                                   as_rvec_array(state_local->x.data()),
6627                                   ncg_moved, bRedist ? comm->movedBuffer.data() : nullptr,
6628                                   fr->nbv->grp[eintLocal].kernel_type,
6629                                   fr->nbv->nbat);
6630
6631                 nbnxn_get_ncells(fr->nbv->nbs.get(), &ncells_new[XX], &ncells_new[YY]);
6632                 break;
6633             case ecutsGROUP:
6634                 fill_grid(&comm->zones, fr->ns->grid, dd->ncg_home,
6635                           0, dd->ncg_home, fr->cg_cm);
6636
6637                 copy_ivec(fr->ns->grid->n, ncells_new);
6638                 break;
6639             default:
6640                 gmx_incons("unimplemented");
6641         }
6642
6643         bResortAll = bMasterState;
6644
6645         /* Check if we can user the old order and ns grid cell indices
6646          * of the charge groups to sort the charge groups efficiently.
6647          */
6648         if (ncells_new[XX] != ncells_old[XX] ||
6649             ncells_new[YY] != ncells_old[YY] ||
6650             ncells_new[ZZ] != ncells_old[ZZ])
6651         {
6652             bResortAll = TRUE;
6653         }
6654
6655         if (debug)
6656         {
6657             fprintf(debug, "Step %s, sorting the %d home charge groups\n",
6658                     gmx_step_str(step, sbuf), dd->ncg_home);
6659         }
6660         dd_sort_state(dd, fr->cg_cm, fr, state_local,
6661                       bResortAll ? -1 : ncg_home_old);
6662
6663         /* After sorting and compacting we set the correct size */
6664         dd_resize_state(state_local, f, comm->atomRanges.numHomeAtoms());
6665
6666         /* Rebuild all the indices */
6667         ga2la_clear(dd->ga2la);
6668         ncgindex_set = 0;
6669
6670         wallcycle_sub_stop(wcycle, ewcsDD_GRID);
6671     }
6672
6673     wallcycle_sub_start(wcycle, ewcsDD_SETUPCOMM);
6674
6675     /* Setup up the communication and communicate the coordinates */
6676     setup_dd_communication(dd, state_local->box, &ddbox, fr, state_local, f);
6677
6678     /* Set the indices */
6679     make_dd_indices(dd, cgs_gl->index, ncgindex_set);
6680
6681     /* Set the charge group boundaries for neighbor searching */
6682     set_cg_boundaries(&comm->zones);
6683
6684     if (fr->cutoff_scheme == ecutsVERLET)
6685     {
6686         /* When bSortCG=true, we have already set the size for zone 0 */
6687         set_zones_size(dd, state_local->box, &ddbox,
6688                        bSortCG ? 1 : 0, comm->zones.n,
6689                        0);
6690     }
6691
6692     wallcycle_sub_stop(wcycle, ewcsDD_SETUPCOMM);
6693
6694     /*
6695        write_dd_pdb("dd_home",step,"dump",top_global,cr,
6696                  -1,as_rvec_array(state_local->x.data()),state_local->box);
6697      */
6698
6699     wallcycle_sub_start(wcycle, ewcsDD_MAKETOP);
6700
6701     /* Extract a local topology from the global topology */
6702     for (int i = 0; i < dd->ndim; i++)
6703     {
6704         np[dd->dim[i]] = comm->cd[i].numPulses();
6705     }
6706     dd_make_local_top(dd, &comm->zones, dd->npbcdim, state_local->box,
6707                       comm->cellsize_min, np,
6708                       fr,
6709                       fr->cutoff_scheme == ecutsGROUP ? fr->cg_cm : as_rvec_array(state_local->x.data()),
6710                       vsite, top_global, top_local);
6711
6712     wallcycle_sub_stop(wcycle, ewcsDD_MAKETOP);
6713
6714     wallcycle_sub_start(wcycle, ewcsDD_MAKECONSTR);
6715
6716     /* Set up the special atom communication */
6717     int n = comm->atomRanges.end(DDAtomRanges::Type::Zones);
6718     for (int i = static_cast<int>(DDAtomRanges::Type::Zones) + 1; i < static_cast<int>(DDAtomRanges::Type::Number); i++)
6719     {
6720         auto range = static_cast<DDAtomRanges::Type>(i);
6721         switch (range)
6722         {
6723             case DDAtomRanges::Type::Vsites:
6724                 if (vsite && vsite->n_intercg_vsite)
6725                 {
6726                     n = dd_make_local_vsites(dd, n, top_local->idef.il);
6727                 }
6728                 break;
6729             case DDAtomRanges::Type::Constraints:
6730                 if (dd->bInterCGcons || dd->bInterCGsettles)
6731                 {
6732                     /* Only for inter-cg constraints we need special code */
6733                     n = dd_make_local_constraints(dd, n, top_global, fr->cginfo,
6734                                                   constr, ir->nProjOrder,
6735                                                   top_local->idef.il);
6736                 }
6737                 break;
6738             default:
6739                 gmx_incons("Unknown special atom type setup");
6740         }
6741         comm->atomRanges.setEnd(range, n);
6742     }
6743
6744     wallcycle_sub_stop(wcycle, ewcsDD_MAKECONSTR);
6745
6746     wallcycle_sub_start(wcycle, ewcsDD_TOPOTHER);
6747
6748     /* Make space for the extra coordinates for virtual site
6749      * or constraint communication.
6750      */
6751     state_local->natoms = comm->atomRanges.numAtomsTotal();
6752
6753     dd_resize_state(state_local, f, state_local->natoms);
6754
6755     if (fr->haveDirectVirialContributions)
6756     {
6757         if (vsite && vsite->n_intercg_vsite)
6758         {
6759             nat_f_novirsum = comm->atomRanges.end(DDAtomRanges::Type::Vsites);
6760         }
6761         else
6762         {
6763             if (EEL_FULL(ir->coulombtype) && dd->n_intercg_excl > 0)
6764             {
6765                 nat_f_novirsum = comm->atomRanges.end(DDAtomRanges::Type::Zones);
6766             }
6767             else
6768             {
6769                 nat_f_novirsum = comm->atomRanges.numHomeAtoms();
6770             }
6771         }
6772     }
6773     else
6774     {
6775         nat_f_novirsum = 0;
6776     }
6777
6778     /* Set the number of atoms required for the force calculation.
6779      * Forces need to be constrained when doing energy
6780      * minimization. For simple simulations we could avoid some
6781      * allocation, zeroing and copying, but this is probably not worth
6782      * the complications and checking.
6783      */
6784     forcerec_set_ranges(fr, dd->ncg_home, dd->globalAtomGroupIndices.size(),
6785                         comm->atomRanges.end(DDAtomRanges::Type::Zones),
6786                         comm->atomRanges.end(DDAtomRanges::Type::Constraints),
6787                         nat_f_novirsum);
6788
6789     /* Update atom data for mdatoms and several algorithms */
6790     mdAlgorithmsSetupAtomData(cr, ir, top_global, top_local, fr,
6791                               nullptr, mdAtoms, constr, vsite, nullptr);
6792
6793     auto mdatoms = mdAtoms->mdatoms();
6794     if (!thisRankHasDuty(cr, DUTY_PME))
6795     {
6796         /* Send the charges and/or c6/sigmas to our PME only node */
6797         gmx_pme_send_parameters(cr,
6798                                 fr->ic,
6799                                 mdatoms->nChargePerturbed, mdatoms->nTypePerturbed,
6800                                 mdatoms->chargeA, mdatoms->chargeB,
6801                                 mdatoms->sqrt_c6A, mdatoms->sqrt_c6B,
6802                                 mdatoms->sigmaA, mdatoms->sigmaB,
6803                                 dd_pme_maxshift_x(dd), dd_pme_maxshift_y(dd));
6804     }
6805
6806     if (ir->bPull)
6807     {
6808         /* Update the local pull groups */
6809         dd_make_local_pull_groups(cr, ir->pull_work);
6810     }
6811
6812     if (ir->bRot)
6813     {
6814         /* Update the local rotation groups */
6815         dd_make_local_rotation_groups(dd, enforcedRotation);
6816     }
6817
6818     if (dd->atomSets != nullptr)
6819     {
6820         /* Update the local atom sets */
6821         dd->atomSets->setIndicesInDomainDecomposition(*(dd->ga2la));
6822     }
6823
6824     /* Update the local atoms to be communicated via the IMD protocol if bIMD is TRUE. */
6825     dd_make_local_IMD_atoms(ir->bIMD, dd, ir->imd);
6826
6827     add_dd_statistics(dd);
6828
6829     /* Make sure we only count the cycles for this DD partitioning */
6830     clear_dd_cycle_counts(dd);
6831
6832     /* Because the order of the atoms might have changed since
6833      * the last vsite construction, we need to communicate the constructing
6834      * atom coordinates again (for spreading the forces this MD step).
6835      */
6836     dd_move_x_vsites(dd, state_local->box, as_rvec_array(state_local->x.data()));
6837
6838     wallcycle_sub_stop(wcycle, ewcsDD_TOPOTHER);
6839
6840     if (comm->nstDDDump > 0 && step % comm->nstDDDump == 0)
6841     {
6842         dd_move_x(dd, state_local->box, state_local->x, nullWallcycle);
6843         write_dd_pdb("dd_dump", step, "dump", top_global, cr,
6844                      -1, as_rvec_array(state_local->x.data()), state_local->box);
6845     }
6846
6847     /* Store the partitioning step */
6848     comm->partition_step = step;
6849
6850     /* Increase the DD partitioning counter */
6851     dd->ddp_count++;
6852     /* The state currently matches this DD partitioning count, store it */
6853     state_local->ddp_count = dd->ddp_count;
6854     if (bMasterState)
6855     {
6856         /* The DD master node knows the complete cg distribution,
6857          * store the count so we can possibly skip the cg info communication.
6858          */
6859         comm->master_cg_ddp_count = (bSortCG ? 0 : dd->ddp_count);
6860     }
6861
6862     if (comm->DD_debug > 0)
6863     {
6864         /* Set the env var GMX_DD_DEBUG if you suspect corrupted indices */
6865         check_index_consistency(dd, top_global->natoms, ncg_mtop(top_global),
6866                                 "after partitioning");
6867     }
6868
6869     wallcycle_stop(wcycle, ewcDOMDEC);
6870 }
6871
6872 /*! \brief Check whether bonded interactions are missing, if appropriate */
6873 void checkNumberOfBondedInteractions(FILE                 *fplog,
6874                                      t_commrec            *cr,
6875                                      int                   totalNumberOfBondedInteractions,
6876                                      const gmx_mtop_t     *top_global,
6877                                      const gmx_localtop_t *top_local,
6878                                      const t_state        *state,
6879                                      bool                 *shouldCheckNumberOfBondedInteractions)
6880 {
6881     if (*shouldCheckNumberOfBondedInteractions)
6882     {
6883         if (totalNumberOfBondedInteractions != cr->dd->nbonded_global)
6884         {
6885             dd_print_missing_interactions(fplog, cr, totalNumberOfBondedInteractions, top_global, top_local, state); // Does not return
6886         }
6887         *shouldCheckNumberOfBondedInteractions = false;
6888     }
6889 }