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