src/gromacs/domdec/domdec.cpp

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