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