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,2019, 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 <cassert>
  43 #include <cinttypes>
  44 #include <climits>
  45 #include <cmath>
  46 #include <cstdio>
  47 #include <cstdlib>
  48 #include <cstring>
  49
  50 #include <algorithm>
  51 #include <memory>
  52
  53 #include "gromacs/domdec/collect.h"
  54 #include "gromacs/domdec/dlb.h"
  55 #include "gromacs/domdec/dlbtiming.h"
  56 #include "gromacs/domdec/domdec_network.h"
  57 #include "gromacs/domdec/ga2la.h"
  58 #include "gromacs/domdec/options.h"
  59 #include "gromacs/domdec/partition.h"
  60 #include "gromacs/gmxlib/network.h"
  61 #include "gromacs/gmxlib/nrnb.h"
  62 #include "gromacs/gpu_utils/gpu_utils.h"
  63 #include "gromacs/hardware/hw_info.h"
  64 #include "gromacs/listed_forces/manage_threading.h"
  65 #include "gromacs/math/vec.h"
  66 #include "gromacs/math/vectypes.h"
  67 #include "gromacs/mdlib/calc_verletbuf.h"
  68 #include "gromacs/mdlib/constr.h"
  69 #include "gromacs/mdlib/constraintrange.h"
  70 #include "gromacs/mdlib/updategroups.h"
  71 #include "gromacs/mdlib/vsite.h"
  72 #include "gromacs/mdtypes/commrec.h"
  73 #include "gromacs/mdtypes/forceoutput.h"
  74 #include "gromacs/mdtypes/inputrec.h"
  75 #include "gromacs/mdtypes/mdrunoptions.h"
  76 #include "gromacs/mdtypes/state.h"
  77 #include "gromacs/pbcutil/ishift.h"
  78 #include "gromacs/pbcutil/pbc.h"
  79 #include "gromacs/pulling/pull.h"
  80 #include "gromacs/timing/wallcycle.h"
  81 #include "gromacs/topology/block.h"
  82 #include "gromacs/topology/idef.h"
  83 #include "gromacs/topology/ifunc.h"
  84 #include "gromacs/topology/mtop_lookup.h"
  85 #include "gromacs/topology/mtop_util.h"
  86 #include "gromacs/topology/topology.h"
  87 #include "gromacs/utility/basedefinitions.h"
  88 #include "gromacs/utility/basenetwork.h"
  89 #include "gromacs/utility/cstringutil.h"
  90 #include "gromacs/utility/exceptions.h"
  91 #include "gromacs/utility/fatalerror.h"
  92 #include "gromacs/utility/gmxmpi.h"
  93 #include "gromacs/utility/logger.h"
  94 #include "gromacs/utility/real.h"
  95 #include "gromacs/utility/smalloc.h"
  96 #include "gromacs/utility/strconvert.h"
  97 #include "gromacs/utility/stringstream.h"
  98 #include "gromacs/utility/stringutil.h"
  99 #include "gromacs/utility/textwriter.h"
 100
 101 #include "atomdistribution.h"
 102 #include "box.h"
 103 #include "cellsizes.h"
 104 #include "distribute.h"
 105 #include "domdec_constraints.h"
 106 #include "domdec_internal.h"
 107 #include "domdec_vsite.h"
 108 #include "redistribute.h"
 109 #include "utility.h"
 110
 111 // TODO remove this when moving domdec into gmx namespace
 112 using gmx::DdRankOrder;
 113 using gmx::DlbOption;
 114 using gmx::DomdecOptions;
 115
 116 static const char *edlbs_names[int(DlbState::nr)] = { "off", "auto", "locked", "on", "on" };
 117
 118 /* The size per atom group of the cggl_flag buffer in gmx_domdec_comm_t */
 119 #define DD_CGIBS 2
 120
 121 /* The flags for the cggl_flag buffer in gmx_domdec_comm_t */
 122 #define DD_FLAG_NRCG  65535
 123 #define DD_FLAG_FW(d) (1<<(16+(d)*2))
 124 #define DD_FLAG_BW(d) (1<<(16+(d)*2+1))
 125
 126 /* The DD zone order */
 127 static const ivec dd_zo[DD_MAXZONE] =
 128 {{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}, {0, 1, 1}, {0, 0, 1}, {1, 0, 1}, {1, 1, 1}};
 129
 130 /* The non-bonded zone-pair setup for domain decomposition
 131  * The first number is the i-zone, the second number the first j-zone seen by
 132  * this i-zone, the third number the last+1 j-zone seen by this i-zone.
 133  * As is, this is for 3D decomposition, where there are 4 i-zones.
 134  * With 2D decomposition use only the first 2 i-zones and a last+1 j-zone of 4.
 135  * With 1D decomposition use only the first i-zone and a last+1 j-zone of 2.
 136  */
 137 static const int
 138     ddNonbondedZonePairRanges[DD_MAXIZONE][3] = {{0, 0, 8},
 139                                                  {1, 3, 6},
 140                                                  {2, 5, 6},
 141                                                  {3, 5, 7}};
 142
 143
 144 /*
 145    #define dd_index(n,i) ((((i)[ZZ]*(n)[YY] + (i)[YY])*(n)[XX]) + (i)[XX])
 146
 147    static void index2xyz(ivec nc,int ind,ivec xyz)
 148    {
 149    xyz[XX] = ind % nc[XX];
 150    xyz[YY] = (ind / nc[XX]) % nc[YY];
 151    xyz[ZZ] = ind / (nc[YY]*nc[XX]);
 152    }
 153  */
 154
 155 static void ddindex2xyz(const ivec nc, int ind, ivec xyz)
 156 {
 157     xyz[XX] = ind / (nc[YY]*nc[ZZ]);
 158     xyz[YY] = (ind / nc[ZZ]) % nc[YY];
 159     xyz[ZZ] = ind % nc[ZZ];
 160 }
 161
 162 static int ddcoord2ddnodeid(gmx_domdec_t *dd, ivec c)
 163 {
 164     int ddindex;
 165     int ddnodeid = -1;
 166
 167     ddindex = dd_index(dd->nc, c);
 168     if (dd->comm->bCartesianPP_PME)
 169     {
 170         ddnodeid = dd->comm->ddindex2ddnodeid[ddindex];
 171     }
 172     else if (dd->comm->bCartesianPP)
 173     {
 174 #if GMX_MPI
 175         MPI_Cart_rank(dd->mpi_comm_all, c, &ddnodeid);
 176 #endif
 177     }
 178     else
 179     {
 180         ddnodeid = ddindex;
 181     }
 182
 183     return ddnodeid;
 184 }
 185
 186 int ddglatnr(const gmx_domdec_t *dd, int i)
 187 {
 188     int atnr;
 189
 190     if (dd == nullptr)
 191     {
 192         atnr = i + 1;
 193     }
 194     else
 195     {
 196         if (i >= dd->comm->atomRanges.numAtomsTotal())
 197         {
 198             gmx_fatal(FARGS, "glatnr called with %d, which is larger than the local number of atoms (%d)", i, dd->comm->atomRanges.numAtomsTotal());
 199         }
 200         atnr = dd->globalAtomIndices[i] + 1;
 201     }
 202
 203     return atnr;
 204 }
 205
 206 t_block *dd_charge_groups_global(gmx_domdec_t *dd)
 207 {
 208     return &dd->comm->cgs_gl;
 209 }
 210
 211 gmx::ArrayRef<const gmx::RangePartitioning> getUpdateGroupingPerMoleculetype(const gmx_domdec_t &dd)
 212 {
 213     GMX_RELEASE_ASSERT(dd.comm, "Need a valid dd.comm");
 214     return dd.comm->systemInfo.updateGroupingPerMoleculetype;
 215 }
 216
 217 void dd_store_state(gmx_domdec_t *dd, t_state *state)
 218 {
 219     int i;
 220
 221     if (state->ddp_count != dd->ddp_count)
 222     {
 223         gmx_incons("The MD state does not match the domain decomposition state");
 224     }
 225
 226     state->cg_gl.resize(dd->ncg_home);
 227     for (i = 0; i < dd->ncg_home; i++)
 228     {
 229         state->cg_gl[i] = dd->globalAtomGroupIndices[i];
 230     }
 231
 232     state->ddp_count_cg_gl = dd->ddp_count;
 233 }
 234
 235 gmx_domdec_zones_t *domdec_zones(gmx_domdec_t *dd)
 236 {
 237     return &dd->comm->zones;
 238 }
 239
 240 void dd_get_ns_ranges(const gmx_domdec_t *dd, int icg,
 241                       int *jcg0, int *jcg1, ivec shift0, ivec shift1)
 242 {
 243     gmx_domdec_zones_t *zones;
 244     int                 izone, d, dim;
 245
 246     zones = &dd->comm->zones;
 247
 248     izone = 0;
 249     while (icg >= zones->izone[izone].cg1)
 250     {
 251         izone++;
 252     }
 253
 254     if (izone == 0)
 255     {
 256         *jcg0 = icg;
 257     }
 258     else if (izone < zones->nizone)
 259     {
 260         *jcg0 = zones->izone[izone].jcg0;
 261     }
 262     else
 263     {
 264         gmx_fatal(FARGS, "DD icg %d out of range: izone (%d) >= nizone (%d)",
 265                   icg, izone, zones->nizone);
 266     }
 267
 268     *jcg1 = zones->izone[izone].jcg1;
 269
 270     for (d = 0; d < dd->ndim; d++)
 271     {
 272         dim         = dd->dim[d];
 273         shift0[dim] = zones->izone[izone].shift0[dim];
 274         shift1[dim] = zones->izone[izone].shift1[dim];
 275         if (dd->comm->tric_dir[dim] || (isDlbOn(dd->comm) && d > 0))
 276         {
 277             /* A conservative approach, this can be optimized */
 278             shift0[dim] -= 1;
 279             shift1[dim] += 1;
 280         }
 281     }
 282 }
 283
 284 int dd_numHomeAtoms(const gmx_domdec_t &dd)
 285 {
 286     return dd.comm->atomRanges.numHomeAtoms();
 287 }
 288
 289 int dd_natoms_mdatoms(const gmx_domdec_t *dd)
 290 {
 291     /* We currently set mdatoms entries for all atoms:
 292      * local + non-local + communicated for vsite + constraints
 293      */
 294
 295     return dd->comm->atomRanges.numAtomsTotal();
 296 }
 297
 298 int dd_natoms_vsite(const gmx_domdec_t *dd)
 299 {
 300     return dd->comm->atomRanges.end(DDAtomRanges::Type::Vsites);
 301 }
 302
 303 void dd_get_constraint_range(const gmx_domdec_t *dd, int *at_start, int *at_end)
 304 {
 305     *at_start = dd->comm->atomRanges.start(DDAtomRanges::Type::Constraints);
 306     *at_end   = dd->comm->atomRanges.end(DDAtomRanges::Type::Constraints);
 307 }
 308
 309 void dd_move_x(gmx_domdec_t             *dd,
 310                const matrix              box,
 311                gmx::ArrayRef<gmx::RVec>  x,
 312                gmx_wallcycle            *wcycle)
 313 {
 314     wallcycle_start(wcycle, ewcMOVEX);
 315
 316     int                    nzone, nat_tot;
 317     gmx_domdec_comm_t     *comm;
 318     gmx_domdec_comm_dim_t *cd;
 319     rvec                   shift = {0, 0, 0};
 320     gmx_bool               bPBC, bScrew;
 321
 322     comm = dd->comm;
 323
 324     nzone   = 1;
 325     nat_tot = comm->atomRanges.numHomeAtoms();
 326     for (int d = 0; d < dd->ndim; d++)
 327     {
 328         bPBC   = (dd->ci[dd->dim[d]] == 0);
 329         bScrew = (bPBC && dd->unitCellInfo.haveScrewPBC && dd->dim[d] == XX);
 330         if (bPBC)
 331         {
 332             copy_rvec(box[dd->dim[d]], shift);
 333         }
 334         cd = &comm->cd[d];
 335         for (const gmx_domdec_ind_t &ind : cd->ind)
 336         {
 337             DDBufferAccess<gmx::RVec>  sendBufferAccess(comm->rvecBuffer, ind.nsend[nzone + 1]);
 338             gmx::ArrayRef<gmx::RVec>  &sendBuffer = sendBufferAccess.buffer;
 339             int                        n          = 0;
 340             if (!bPBC)
 341             {
 342                 for (int j : ind.index)
 343                 {
 344                     sendBuffer[n] = x[j];
 345                     n++;
 346                 }
 347             }
 348             else if (!bScrew)
 349             {
 350                 for (int j : ind.index)
 351                 {
 352                     /* We need to shift the coordinates */
 353                     for (int d = 0; d < DIM; d++)
 354                     {
 355                         sendBuffer[n][d] = x[j][d] + shift[d];
 356                     }
 357                     n++;
 358                 }
 359             }
 360             else
 361             {
 362                 for (int j : ind.index)
 363                 {
 364                     /* Shift x */
 365                     sendBuffer[n][XX] = x[j][XX] + shift[XX];
 366                     /* Rotate y and z.
 367                      * This operation requires a special shift force
 368                      * treatment, which is performed in calc_vir.
 369                      */
 370                     sendBuffer[n][YY] = box[YY][YY] - x[j][YY];
 371                     sendBuffer[n][ZZ] = box[ZZ][ZZ] - x[j][ZZ];
 372                     n++;
 373                 }
 374             }
 375
 376             DDBufferAccess<gmx::RVec>  receiveBufferAccess(comm->rvecBuffer2, cd->receiveInPlace ? 0 : ind.nrecv[nzone + 1]);
 377
 378             gmx::ArrayRef<gmx::RVec>   receiveBuffer;
 379             if (cd->receiveInPlace)
 380             {
 381                 receiveBuffer = gmx::arrayRefFromArray(x.data() + nat_tot, ind.nrecv[nzone + 1]);
 382             }
 383             else
 384             {
 385                 receiveBuffer = receiveBufferAccess.buffer;
 386             }
 387             /* Send and receive the coordinates */
 388             ddSendrecv(dd, d, dddirBackward,
 389                        sendBuffer, receiveBuffer);
 390
 391             if (!cd->receiveInPlace)
 392             {
 393                 int j = 0;
 394                 for (int zone = 0; zone < nzone; zone++)
 395                 {
 396                     for (int i = ind.cell2at0[zone]; i < ind.cell2at1[zone]; i++)
 397                     {
 398                         x[i] = receiveBuffer[j++];
 399                     }
 400                 }
 401             }
 402             nat_tot += ind.nrecv[nzone+1];
 403         }
 404         nzone += nzone;
 405     }
 406
 407     wallcycle_stop(wcycle, ewcMOVEX);
 408 }
 409
 410 void dd_move_f(gmx_domdec_t              *dd,
 411                gmx::ForceWithShiftForces *forceWithShiftForces,
 412                gmx_wallcycle             *wcycle)
 413 {
 414     wallcycle_start(wcycle, ewcMOVEF);
 415
 416     gmx::ArrayRef<gmx::RVec> f       = forceWithShiftForces->force();
 417     gmx::ArrayRef<gmx::RVec> fshift  = forceWithShiftForces->shiftForces();
 418
 419     gmx_domdec_comm_t       &comm    = *dd->comm;
 420     int                      nzone   = comm.zones.n/2;
 421     int                      nat_tot = comm.atomRanges.end(DDAtomRanges::Type::Zones);
 422     for (int d = dd->ndim-1; d >= 0; d--)
 423     {
 424         /* Only forces in domains near the PBC boundaries need to
 425            consider PBC in the treatment of fshift */
 426         const bool shiftForcesNeedPbc = (forceWithShiftForces->computeVirial() && dd->ci[dd->dim[d]] == 0);
 427         const bool applyScrewPbc      = (shiftForcesNeedPbc && dd->unitCellInfo.haveScrewPBC && dd->dim[d] == XX);
 428         /* Determine which shift vector we need */
 429         ivec       vis   = { 0, 0, 0 };
 430         vis[dd->dim[d]]  = 1;
 431         const int  is    = IVEC2IS(vis);
 432
 433         /* Loop over the pulses */
 434         const gmx_domdec_comm_dim_t &cd = comm.cd[d];
 435         for (int p = cd.numPulses() - 1; p >= 0; p--)
 436         {
 437             const gmx_domdec_ind_t    &ind  = cd.ind[p];
 438             DDBufferAccess<gmx::RVec>  receiveBufferAccess(comm.rvecBuffer, ind.nsend[nzone + 1]);
 439             gmx::ArrayRef<gmx::RVec>  &receiveBuffer = receiveBufferAccess.buffer;
 440
 441             nat_tot                                 -= ind.nrecv[nzone+1];
 442
 443             DDBufferAccess<gmx::RVec>  sendBufferAccess(comm.rvecBuffer2, cd.receiveInPlace ? 0 : ind.nrecv[nzone + 1]);
 444
 445             gmx::ArrayRef<gmx::RVec>   sendBuffer;
 446             if (cd.receiveInPlace)
 447             {
 448                 sendBuffer = gmx::arrayRefFromArray(f.data() + nat_tot, ind.nrecv[nzone + 1]);
 449             }
 450             else
 451             {
 452                 sendBuffer = sendBufferAccess.buffer;
 453                 int j = 0;
 454                 for (int zone = 0; zone < nzone; zone++)
 455                 {
 456                     for (int i = ind.cell2at0[zone]; i < ind.cell2at1[zone]; i++)
 457                     {
 458                         sendBuffer[j++] = f[i];
 459                     }
 460                 }
 461             }
 462             /* Communicate the forces */
 463             ddSendrecv(dd, d, dddirForward,
 464                        sendBuffer, receiveBuffer);
 465             /* Add the received forces */
 466             int n = 0;
 467             if (!shiftForcesNeedPbc)
 468             {
 469                 for (int j : ind.index)
 470                 {
 471                     for (int d = 0; d < DIM; d++)
 472                     {
 473                         f[j][d] += receiveBuffer[n][d];
 474                     }
 475                     n++;
 476                 }
 477             }
 478             else if (!applyScrewPbc)
 479             {
 480                 for (int j : ind.index)
 481                 {
 482                     for (int d = 0; d < DIM; d++)
 483                     {
 484                         f[j][d] += receiveBuffer[n][d];
 485                     }
 486                     /* Add this force to the shift force */
 487                     for (int d = 0; d < DIM; d++)
 488                     {
 489                         fshift[is][d] += receiveBuffer[n][d];
 490                     }
 491                     n++;
 492                 }
 493             }
 494             else
 495             {
 496                 for (int j : ind.index)
 497                 {
 498                     /* Rotate the force */
 499                     f[j][XX] += receiveBuffer[n][XX];
 500                     f[j][YY] -= receiveBuffer[n][YY];
 501                     f[j][ZZ] -= receiveBuffer[n][ZZ];
 502                     if (shiftForcesNeedPbc)
 503                     {
 504                         /* Add this force to the shift force */
 505                         for (int d = 0; d < DIM; d++)
 506                         {
 507                             fshift[is][d] += receiveBuffer[n][d];
 508                         }
 509                     }
 510                     n++;
 511                 }
 512             }
 513         }
 514         nzone /= 2;
 515     }
 516     wallcycle_stop(wcycle, ewcMOVEF);
 517 }
 518
 519 /* Convenience function for extracting a real buffer from an rvec buffer
 520  *
 521  * To reduce the number of temporary communication buffers and avoid
 522  * cache polution, we reuse gmx::RVec buffers for storing reals.
 523  * This functions return a real buffer reference with the same number
 524  * of elements as the gmx::RVec buffer (so 1/3 of the size in bytes).
 525  */
 526 static gmx::ArrayRef<real>
 527 realArrayRefFromRvecArrayRef(gmx::ArrayRef<gmx::RVec> arrayRef)
 528 {
 529     return gmx::arrayRefFromArray(as_rvec_array(arrayRef.data())[0],
 530                                   arrayRef.size());
 531 }
 532
 533 void dd_atom_spread_real(gmx_domdec_t *dd, real v[])
 534 {
 535     int                    nzone, nat_tot;
 536     gmx_domdec_comm_t     *comm;
 537     gmx_domdec_comm_dim_t *cd;
 538
 539     comm = dd->comm;
 540
 541     nzone   = 1;
 542     nat_tot = comm->atomRanges.numHomeAtoms();
 543     for (int d = 0; d < dd->ndim; d++)
 544     {
 545         cd = &comm->cd[d];
 546         for (const gmx_domdec_ind_t &ind : cd->ind)
 547         {
 548             /* Note: We provision for RVec instead of real, so a factor of 3
 549              * more than needed. The buffer actually already has this size
 550              * and we pass a plain pointer below, so this does not matter.
 551              */
 552             DDBufferAccess<gmx::RVec> sendBufferAccess(comm->rvecBuffer, ind.nsend[nzone + 1]);
 553             gmx::ArrayRef<real>       sendBuffer = realArrayRefFromRvecArrayRef(sendBufferAccess.buffer);
 554             int                       n          = 0;
 555             for (int j : ind.index)
 556             {
 557                 sendBuffer[n++] = v[j];
 558             }
 559
 560             DDBufferAccess<gmx::RVec> receiveBufferAccess(comm->rvecBuffer2, cd->receiveInPlace ? 0 : ind.nrecv[nzone + 1]);
 561
 562             gmx::ArrayRef<real>       receiveBuffer;
 563             if (cd->receiveInPlace)
 564             {
 565                 receiveBuffer = gmx::arrayRefFromArray(v + nat_tot, ind.nrecv[nzone + 1]);
 566             }
 567             else
 568             {
 569                 receiveBuffer = realArrayRefFromRvecArrayRef(receiveBufferAccess.buffer);
 570             }
 571             /* Send and receive the data */
 572             ddSendrecv(dd, d, dddirBackward,
 573                        sendBuffer, receiveBuffer);
 574             if (!cd->receiveInPlace)
 575             {
 576                 int j = 0;
 577                 for (int zone = 0; zone < nzone; zone++)
 578                 {
 579                     for (int i = ind.cell2at0[zone]; i < ind.cell2at1[zone]; i++)
 580                     {
 581                         v[i] = receiveBuffer[j++];
 582                     }
 583                 }
 584             }
 585             nat_tot += ind.nrecv[nzone+1];
 586         }
 587         nzone += nzone;
 588     }
 589 }
 590
 591 void dd_atom_sum_real(gmx_domdec_t *dd, real v[])
 592 {
 593     int                    nzone, nat_tot;
 594     gmx_domdec_comm_t     *comm;
 595     gmx_domdec_comm_dim_t *cd;
 596
 597     comm = dd->comm;
 598
 599     nzone   = comm->zones.n/2;
 600     nat_tot = comm->atomRanges.end(DDAtomRanges::Type::Zones);
 601     for (int d = dd->ndim-1; d >= 0; d--)
 602     {
 603         cd = &comm->cd[d];
 604         for (int p = cd->numPulses() - 1; p >= 0; p--)
 605         {
 606             const gmx_domdec_ind_t &ind = cd->ind[p];
 607
 608             /* Note: We provision for RVec instead of real, so a factor of 3
 609              * more than needed. The buffer actually already has this size
 610              * and we typecast, so this works as intended.
 611              */
 612             DDBufferAccess<gmx::RVec> receiveBufferAccess(comm->rvecBuffer, ind.nsend[nzone + 1]);
 613             gmx::ArrayRef<real>       receiveBuffer = realArrayRefFromRvecArrayRef(receiveBufferAccess.buffer);
 614             nat_tot -= ind.nrecv[nzone + 1];
 615
 616             DDBufferAccess<gmx::RVec> sendBufferAccess(comm->rvecBuffer2, cd->receiveInPlace ? 0 : ind.nrecv[nzone + 1]);
 617
 618             gmx::ArrayRef<real>       sendBuffer;
 619             if (cd->receiveInPlace)
 620             {
 621                 sendBuffer = gmx::arrayRefFromArray(v + nat_tot, ind.nrecv[nzone + 1]);
 622             }
 623             else
 624             {
 625                 sendBuffer = realArrayRefFromRvecArrayRef(sendBufferAccess.buffer);
 626                 int j = 0;
 627                 for (int zone = 0; zone < nzone; zone++)
 628                 {
 629                     for (int i = ind.cell2at0[zone]; i < ind.cell2at1[zone]; i++)
 630                     {
 631                         sendBuffer[j++] = v[i];
 632                     }
 633                 }
 634             }
 635             /* Communicate the forces */
 636             ddSendrecv(dd, d, dddirForward,
 637                        sendBuffer, receiveBuffer);
 638             /* Add the received forces */
 639             int n = 0;
 640             for (int j : ind.index)
 641             {
 642                 v[j] += receiveBuffer[n];
 643                 n++;
 644             }
 645         }
 646         nzone /= 2;
 647     }
 648 }
 649
 650 real dd_cutoff_multibody(const gmx_domdec_t *dd)
 651 {
 652     const gmx_domdec_comm_t &comm       = *dd->comm;
 653     const DDSystemInfo      &systemInfo = comm.systemInfo;
 654
 655     real                     r = -1;
 656     if (systemInfo.haveInterDomainMultiBodyBondeds)
 657     {
 658         if (comm.cutoff_mbody > 0)
 659         {
 660             r = comm.cutoff_mbody;
 661         }
 662         else
 663         {
 664             /* cutoff_mbody=0 means we do not have DLB */
 665             r = comm.cellsize_min[dd->dim[0]];
 666             for (int di = 1; di < dd->ndim; di++)
 667             {
 668                 r = std::min(r, comm.cellsize_min[dd->dim[di]]);
 669             }
 670             if (comm.systemInfo.filterBondedCommunication)
 671             {
 672                 r = std::max(r, comm.cutoff_mbody);
 673             }
 674             else
 675             {
 676                 r = std::min(r, systemInfo.cutoff);
 677             }
 678         }
 679     }
 680
 681     return r;
 682 }
 683
 684 real dd_cutoff_twobody(const gmx_domdec_t *dd)
 685 {
 686     real r_mb;
 687
 688     r_mb = dd_cutoff_multibody(dd);
 689
 690     return std::max(dd->comm->systemInfo.cutoff, r_mb);
 691 }
 692
 693
 694 static void dd_cart_coord2pmecoord(const gmx_domdec_t *dd, const ivec coord,
 695                                    ivec coord_pme)
 696 {
 697     int nc, ntot;
 698
 699     nc   = dd->nc[dd->comm->cartpmedim];
 700     ntot = dd->comm->ntot[dd->comm->cartpmedim];
 701     copy_ivec(coord, coord_pme);
 702     coord_pme[dd->comm->cartpmedim] =
 703         nc + (coord[dd->comm->cartpmedim]*(ntot - nc) + (ntot - nc)/2)/nc;
 704 }
 705
 706 static int ddindex2pmeindex(const gmx_domdec_t *dd, int ddindex)
 707 {
 708     int npp, npme;
 709
 710     npp  = dd->nnodes;
 711     npme = dd->comm->npmenodes;
 712
 713     /* Here we assign a PME node to communicate with this DD node
 714      * by assuming that the major index of both is x.
 715      * We add cr->npmenodes/2 to obtain an even distribution.
 716      */
 717     return (ddindex*npme + npme/2)/npp;
 718 }
 719
 720 static int *dd_interleaved_pme_ranks(const gmx_domdec_t *dd)
 721 {
 722     int *pme_rank;
 723     int  n, i, p0, p1;
 724
 725     snew(pme_rank, dd->comm->npmenodes);
 726     n = 0;
 727     for (i = 0; i < dd->nnodes; i++)
 728     {
 729         p0 = ddindex2pmeindex(dd, i);
 730         p1 = ddindex2pmeindex(dd, i+1);
 731         if (i+1 == dd->nnodes || p1 > p0)
 732         {
 733             if (debug)
 734             {
 735                 fprintf(debug, "pme_rank[%d] = %d\n", n, i+1+n);
 736             }
 737             pme_rank[n] = i + 1 + n;
 738             n++;
 739         }
 740     }
 741
 742     return pme_rank;
 743 }
 744
 745 static int gmx_ddcoord2pmeindex(const t_commrec *cr, int x, int y, int z)
 746 {
 747     gmx_domdec_t *dd;
 748     ivec          coords;
 749     int           slab;
 750
 751     dd = cr->dd;
 752     /*
 753        if (dd->comm->bCartesian) {
 754        gmx_ddindex2xyz(dd->nc,ddindex,coords);
 755        dd_coords2pmecoords(dd,coords,coords_pme);
 756        copy_ivec(dd->ntot,nc);
 757        nc[dd->cartpmedim]         -= dd->nc[dd->cartpmedim];
 758        coords_pme[dd->cartpmedim] -= dd->nc[dd->cartpmedim];
 759
 760        slab = (coords_pme[XX]*nc[YY] + coords_pme[YY])*nc[ZZ] + coords_pme[ZZ];
 761        } else {
 762        slab = (ddindex*cr->npmenodes + cr->npmenodes/2)/dd->nnodes;
 763        }
 764      */
 765     coords[XX] = x;
 766     coords[YY] = y;
 767     coords[ZZ] = z;
 768     slab       = ddindex2pmeindex(dd, dd_index(dd->nc, coords));
 769
 770     return slab;
 771 }
 772
 773 static int ddcoord2simnodeid(const t_commrec *cr, int x, int y, int z)
 774 {
 775     gmx_domdec_comm_t *comm;
 776     ivec               coords;
 777     int                ddindex, nodeid = -1;
 778
 779     comm = cr->dd->comm;
 780
 781     coords[XX] = x;
 782     coords[YY] = y;
 783     coords[ZZ] = z;
 784     if (comm->bCartesianPP_PME)
 785     {
 786 #if GMX_MPI
 787         MPI_Cart_rank(cr->mpi_comm_mysim, coords, &nodeid);
 788 #endif
 789     }
 790     else
 791     {
 792         ddindex = dd_index(cr->dd->nc, coords);
 793         if (comm->bCartesianPP)
 794         {
 795             nodeid = comm->ddindex2simnodeid[ddindex];
 796         }
 797         else
 798         {
 799             if (comm->pmenodes)
 800             {
 801                 nodeid = ddindex + gmx_ddcoord2pmeindex(cr, x, y, z);
 802             }
 803             else
 804             {
 805                 nodeid = ddindex;
 806             }
 807         }
 808     }
 809
 810     return nodeid;
 811 }
 812
 813 static int dd_simnode2pmenode(const gmx_domdec_t         *dd,
 814                               const t_commrec gmx_unused *cr,
 815                               int                         sim_nodeid)
 816 {
 817     int pmenode = -1;
 818
 819     const gmx_domdec_comm_t *comm = dd->comm;
 820
 821     /* This assumes a uniform x domain decomposition grid cell size */
 822     if (comm->bCartesianPP_PME)
 823     {
 824 #if GMX_MPI
 825         ivec coord, coord_pme;
 826         MPI_Cart_coords(cr->mpi_comm_mysim, sim_nodeid, DIM, coord);
 827         if (coord[comm->cartpmedim] < dd->nc[comm->cartpmedim])
 828         {
 829             /* This is a PP node */
 830             dd_cart_coord2pmecoord(dd, coord, coord_pme);
 831             MPI_Cart_rank(cr->mpi_comm_mysim, coord_pme, &pmenode);
 832         }
 833 #endif
 834     }
 835     else if (comm->bCartesianPP)
 836     {
 837         if (sim_nodeid < dd->nnodes)
 838         {
 839             pmenode = dd->nnodes + ddindex2pmeindex(dd, sim_nodeid);
 840         }
 841     }
 842     else
 843     {
 844         /* This assumes DD cells with identical x coordinates
 845          * are numbered sequentially.
 846          */
 847         if (dd->comm->pmenodes == nullptr)
 848         {
 849             if (sim_nodeid < dd->nnodes)
 850             {
 851                 /* The DD index equals the nodeid */
 852                 pmenode = dd->nnodes + ddindex2pmeindex(dd, sim_nodeid);
 853             }
 854         }
 855         else
 856         {
 857             int i = 0;
 858             while (sim_nodeid > dd->comm->pmenodes[i])
 859             {
 860                 i++;
 861             }
 862             if (sim_nodeid < dd->comm->pmenodes[i])
 863             {
 864                 pmenode = dd->comm->pmenodes[i];
 865             }
 866         }
 867     }
 868
 869     return pmenode;
 870 }
 871
 872 NumPmeDomains getNumPmeDomains(const gmx_domdec_t *dd)
 873 {
 874     if (dd != nullptr)
 875     {
 876         return {
 877                    dd->comm->npmenodes_x, dd->comm->npmenodes_y
 878         };
 879     }
 880     else
 881     {
 882         return {
 883                    1, 1
 884         };
 885     }
 886 }
 887
 888 std::vector<int> get_pme_ddranks(const t_commrec *cr, int pmenodeid)
 889 {
 890     gmx_domdec_t *dd;
 891     int           x, y, z;
 892     ivec          coord, coord_pme;
 893
 894     dd = cr->dd;
 895
 896     std::vector<int> ddranks;
 897     ddranks.reserve((dd->nnodes+cr->npmenodes-1)/cr->npmenodes);
 898
 899     for (x = 0; x < dd->nc[XX]; x++)
 900     {
 901         for (y = 0; y < dd->nc[YY]; y++)
 902         {
 903             for (z = 0; z < dd->nc[ZZ]; z++)
 904             {
 905                 if (dd->comm->bCartesianPP_PME)
 906                 {
 907                     coord[XX] = x;
 908                     coord[YY] = y;
 909                     coord[ZZ] = z;
 910                     dd_cart_coord2pmecoord(dd, coord, coord_pme);
 911                     if (dd->ci[XX] == coord_pme[XX] &&
 912                         dd->ci[YY] == coord_pme[YY] &&
 913                         dd->ci[ZZ] == coord_pme[ZZ])
 914                     {
 915                         ddranks.push_back(ddcoord2simnodeid(cr, x, y, z));
 916                     }
 917                 }
 918                 else
 919                 {
 920                     /* The slab corresponds to the nodeid in the PME group */
 921                     if (gmx_ddcoord2pmeindex(cr, x, y, z) == pmenodeid)
 922                     {
 923                         ddranks.push_back(ddcoord2simnodeid(cr, x, y, z));
 924                     }
 925                 }
 926             }
 927         }
 928     }
 929     return ddranks;
 930 }
 931
 932 static gmx_bool receive_vir_ener(const gmx_domdec_t *dd, const t_commrec *cr)
 933 {
 934     gmx_bool bReceive = TRUE;
 935
 936     if (cr->npmenodes < dd->nnodes)
 937     {
 938         gmx_domdec_comm_t *comm = dd->comm;
 939         if (comm->bCartesianPP_PME)
 940         {
 941 #if GMX_MPI
 942             int  pmenode = dd_simnode2pmenode(dd, cr, cr->sim_nodeid);
 943             ivec coords;
 944             MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, coords);
 945             coords[comm->cartpmedim]++;
 946             if (coords[comm->cartpmedim] < dd->nc[comm->cartpmedim])
 947             {
 948                 int rank;
 949                 MPI_Cart_rank(cr->mpi_comm_mysim, coords, &rank);
 950                 if (dd_simnode2pmenode(dd, cr, rank) == pmenode)
 951                 {
 952                     /* This is not the last PP node for pmenode */
 953                     bReceive = FALSE;
 954                 }
 955             }
 956 #else
 957             GMX_RELEASE_ASSERT(false, "Without MPI we should not have Cartesian PP-PME with #PMEnodes < #DDnodes");
 958 #endif
 959         }
 960         else
 961         {
 962             int pmenode = dd_simnode2pmenode(dd, cr, cr->sim_nodeid);
 963             if (cr->sim_nodeid+1 < cr->nnodes &&
 964                 dd_simnode2pmenode(dd, cr, cr->sim_nodeid+1) == pmenode)
 965             {
 966                 /* This is not the last PP node for pmenode */
 967                 bReceive = FALSE;
 968             }
 969         }
 970     }
 971
 972     return bReceive;
 973 }
 974
 975 static void set_slb_pme_dim_f(gmx_domdec_t *dd, int dim, real **dim_f)
 976 {
 977     gmx_domdec_comm_t *comm;
 978     int                i;
 979
 980     comm = dd->comm;
 981
 982     snew(*dim_f, dd->nc[dim]+1);
 983     (*dim_f)[0] = 0;
 984     for (i = 1; i < dd->nc[dim]; i++)
 985     {
 986         if (comm->slb_frac[dim])
 987         {
 988             (*dim_f)[i] = (*dim_f)[i-1] + comm->slb_frac[dim][i-1];
 989         }
 990         else
 991         {
 992             (*dim_f)[i] = static_cast<real>(i)/static_cast<real>(dd->nc[dim]);
 993         }
 994     }
 995     (*dim_f)[dd->nc[dim]] = 1;
 996 }
 997
 998 static void init_ddpme(gmx_domdec_t *dd, gmx_ddpme_t *ddpme, int dimind)
 999 {
1000     int  pmeindex, slab, nso, i;
1001     ivec xyz;
1002
1003     if (dimind == 0 && dd->dim[0] == YY && dd->comm->npmenodes_x == 1)
1004     {
1005         ddpme->dim = YY;
1006     }
1007     else
1008     {
1009         ddpme->dim = dimind;
1010     }
1011     ddpme->dim_match = (ddpme->dim == dd->dim[dimind]);
1012
1013     ddpme->nslab = (ddpme->dim == 0 ?
1014                     dd->comm->npmenodes_x :
1015                     dd->comm->npmenodes_y);
1016
1017     if (ddpme->nslab <= 1)
1018     {
1019         return;
1020     }
1021
1022     nso = dd->comm->npmenodes/ddpme->nslab;
1023     /* Determine for each PME slab the PP location range for dimension dim */
1024     snew(ddpme->pp_min, ddpme->nslab);
1025     snew(ddpme->pp_max, ddpme->nslab);
1026     for (slab = 0; slab < ddpme->nslab; slab++)
1027     {
1028         ddpme->pp_min[slab] = dd->nc[dd->dim[dimind]] - 1;
1029         ddpme->pp_max[slab] = 0;
1030     }
1031     for (i = 0; i < dd->nnodes; i++)
1032     {
1033         ddindex2xyz(dd->nc, i, xyz);
1034         /* For y only use our y/z slab.
1035          * This assumes that the PME x grid size matches the DD grid size.
1036          */
1037         if (dimind == 0 || xyz[XX] == dd->ci[XX])
1038         {
1039             pmeindex = ddindex2pmeindex(dd, i);
1040             if (dimind == 0)
1041             {
1042                 slab = pmeindex/nso;
1043             }
1044             else
1045             {
1046                 slab = pmeindex % ddpme->nslab;
1047             }
1048             ddpme->pp_min[slab] = std::min(ddpme->pp_min[slab], xyz[dimind]);
1049             ddpme->pp_max[slab] = std::max(ddpme->pp_max[slab], xyz[dimind]);
1050         }
1051     }
1052
1053     set_slb_pme_dim_f(dd, ddpme->dim, &ddpme->slb_dim_f);
1054 }
1055
1056 int dd_pme_maxshift_x(const gmx_domdec_t *dd)
1057 {
1058     if (dd->comm->ddpme[0].dim == XX)
1059     {
1060         return dd->comm->ddpme[0].maxshift;
1061     }
1062     else
1063     {
1064         return 0;
1065     }
1066 }
1067
1068 int dd_pme_maxshift_y(const gmx_domdec_t *dd)
1069 {
1070     if (dd->comm->ddpme[0].dim == YY)
1071     {
1072         return dd->comm->ddpme[0].maxshift;
1073     }
1074     else if (dd->comm->npmedecompdim >= 2 && dd->comm->ddpme[1].dim == YY)
1075     {
1076         return dd->comm->ddpme[1].maxshift;
1077     }
1078     else
1079     {
1080         return 0;
1081     }
1082 }
1083
1084 bool ddHaveSplitConstraints(const gmx_domdec_t &dd)
1085 {
1086     return dd.comm->systemInfo.haveSplitConstraints;
1087 }
1088
1089 void dd_cycles_add(const gmx_domdec_t *dd, float cycles, int ddCycl)
1090 {
1091     /* Note that the cycles value can be incorrect, either 0 or some
1092      * extremely large value, when our thread migrated to another core
1093      * with an unsynchronized cycle counter. If this happens less often
1094      * that once per nstlist steps, this will not cause issues, since
1095      * we later subtract the maximum value from the sum over nstlist steps.
1096      * A zero count will slightly lower the total, but that's a small effect.
1097      * Note that the main purpose of the subtraction of the maximum value
1098      * is to avoid throwing off the load balancing when stalls occur due
1099      * e.g. system activity or network congestion.
1100      */
1101     dd->comm->cycl[ddCycl] += cycles;
1102     dd->comm->cycl_n[ddCycl]++;
1103     if (cycles > dd->comm->cycl_max[ddCycl])
1104     {
1105         dd->comm->cycl_max[ddCycl] = cycles;
1106     }
1107 }
1108
1109 #if GMX_MPI
1110 static void make_load_communicator(gmx_domdec_t *dd, int dim_ind, ivec loc)
1111 {
1112     MPI_Comm           c_row;
1113     int                dim, i, rank;
1114     ivec               loc_c;
1115     gmx_bool           bPartOfGroup = FALSE;
1116
1117     dim = dd->dim[dim_ind];
1118     copy_ivec(loc, loc_c);
1119     for (i = 0; i < dd->nc[dim]; i++)
1120     {
1121         loc_c[dim] = i;
1122         rank       = dd_index(dd->nc, loc_c);
1123         if (rank == dd->rank)
1124         {
1125             /* This process is part of the group */
1126             bPartOfGroup = TRUE;
1127         }
1128     }
1129     MPI_Comm_split(dd->mpi_comm_all, bPartOfGroup ? 0 : MPI_UNDEFINED, dd->rank,
1130                    &c_row);
1131     if (bPartOfGroup)
1132     {
1133         dd->comm->mpi_comm_load[dim_ind] = c_row;
1134         if (!isDlbDisabled(dd->comm))
1135         {
1136             DDCellsizesWithDlb &cellsizes = dd->comm->cellsizesWithDlb[dim_ind];
1137
1138             if (dd->ci[dim] == dd->master_ci[dim])
1139             {
1140                 /* This is the root process of this row */
1141                 cellsizes.rowMaster  = std::make_unique<RowMaster>();
1142
1143                 RowMaster &rowMaster = *cellsizes.rowMaster;
1144                 rowMaster.cellFrac.resize(ddCellFractionBufferSize(dd, dim_ind));
1145                 rowMaster.oldCellFrac.resize(dd->nc[dim] + 1);
1146                 rowMaster.isCellMin.resize(dd->nc[dim]);
1147                 if (dim_ind > 0)
1148                 {
1149                     rowMaster.bounds.resize(dd->nc[dim]);
1150                 }
1151                 rowMaster.buf_ncd.resize(dd->nc[dim]);
1152             }
1153             else
1154             {
1155                 /* This is not a root process, we only need to receive cell_f */
1156                 cellsizes.fracRow.resize(ddCellFractionBufferSize(dd, dim_ind));
1157             }
1158         }
1159         if (dd->ci[dim] == dd->master_ci[dim])
1160         {
1161             snew(dd->comm->load[dim_ind].load, dd->nc[dim]*DD_NLOAD_MAX);
1162         }
1163     }
1164 }
1165 #endif
1166
1167 void dd_setup_dlb_resource_sharing(t_commrec            *cr,
1168                                    int                   gpu_id)
1169 {
1170 #if GMX_MPI
1171     int           physicalnode_id_hash;
1172     gmx_domdec_t *dd;
1173     MPI_Comm      mpi_comm_pp_physicalnode;
1174
1175     if (!thisRankHasDuty(cr, DUTY_PP) || gpu_id < 0)
1176     {
1177         /* Only ranks with short-ranged tasks (currently) use GPUs.
1178          * If we don't have GPUs assigned, there are no resources to share.
1179          */
1180         return;
1181     }
1182
1183     physicalnode_id_hash = gmx_physicalnode_id_hash();
1184
1185     dd = cr->dd;
1186
1187     if (debug)
1188     {
1189         fprintf(debug, "dd_setup_dd_dlb_gpu_sharing:\n");
1190         fprintf(debug, "DD PP rank %d physical node hash %d gpu_id %d\n",
1191                 dd->rank, physicalnode_id_hash, gpu_id);
1192     }
1193     /* Split the PP communicator over the physical nodes */
1194     /* TODO: See if we should store this (before), as it's also used for
1195      * for the nodecomm summation.
1196      */
1197     // TODO PhysicalNodeCommunicator could be extended/used to handle
1198     // the need for per-node per-group communicators.
1199     MPI_Comm_split(dd->mpi_comm_all, physicalnode_id_hash, dd->rank,
1200                    &mpi_comm_pp_physicalnode);
1201     MPI_Comm_split(mpi_comm_pp_physicalnode, gpu_id, dd->rank,
1202                    &dd->comm->mpi_comm_gpu_shared);
1203     MPI_Comm_free(&mpi_comm_pp_physicalnode);
1204     MPI_Comm_size(dd->comm->mpi_comm_gpu_shared, &dd->comm->nrank_gpu_shared);
1205
1206     if (debug)
1207     {
1208         fprintf(debug, "nrank_gpu_shared %d\n", dd->comm->nrank_gpu_shared);
1209     }
1210
1211     /* Note that some ranks could share a GPU, while others don't */
1212
1213     if (dd->comm->nrank_gpu_shared == 1)
1214     {
1215         MPI_Comm_free(&dd->comm->mpi_comm_gpu_shared);
1216     }
1217 #else
1218     GMX_UNUSED_VALUE(cr);
1219     GMX_UNUSED_VALUE(gpu_id);
1220 #endif
1221 }
1222
1223 static void make_load_communicators(gmx_domdec_t gmx_unused *dd)
1224 {
1225 #if GMX_MPI
1226     int  dim0, dim1, i, j;
1227     ivec loc;
1228
1229     if (debug)
1230     {
1231         fprintf(debug, "Making load communicators\n");
1232     }
1233
1234     dd->comm->load = new domdec_load_t[std::max(dd->ndim, 1)];
1235     snew(dd->comm->mpi_comm_load, std::max(dd->ndim, 1));
1236
1237     if (dd->ndim == 0)
1238     {
1239         return;
1240     }
1241
1242     clear_ivec(loc);
1243     make_load_communicator(dd, 0, loc);
1244     if (dd->ndim > 1)
1245     {
1246         dim0 = dd->dim[0];
1247         for (i = 0; i < dd->nc[dim0]; i++)
1248         {
1249             loc[dim0] = i;
1250             make_load_communicator(dd, 1, loc);
1251         }
1252     }
1253     if (dd->ndim > 2)
1254     {
1255         dim0 = dd->dim[0];
1256         for (i = 0; i < dd->nc[dim0]; i++)
1257         {
1258             loc[dim0] = i;
1259             dim1      = dd->dim[1];
1260             for (j = 0; j < dd->nc[dim1]; j++)
1261             {
1262                 loc[dim1] = j;
1263                 make_load_communicator(dd, 2, loc);
1264             }
1265         }
1266     }
1267
1268     if (debug)
1269     {
1270         fprintf(debug, "Finished making load communicators\n");
1271     }
1272 #endif
1273 }
1274
1275 /*! \brief Sets up the relation between neighboring domains and zones */
1276 static void setup_neighbor_relations(gmx_domdec_t *dd)
1277 {
1278     int                     d, dim, i, j, m;
1279     ivec                    tmp, s;
1280     gmx_domdec_zones_t     *zones;
1281     gmx_domdec_ns_ranges_t *izone;
1282     GMX_ASSERT(dd->ndim >= 0, "Must have non-negative number of dimensions for DD");
1283
1284     for (d = 0; d < dd->ndim; d++)
1285     {
1286         dim = dd->dim[d];
1287         copy_ivec(dd->ci, tmp);
1288         tmp[dim]           = (tmp[dim] + 1) % dd->nc[dim];
1289         dd->neighbor[d][0] = ddcoord2ddnodeid(dd, tmp);
1290         copy_ivec(dd->ci, tmp);
1291         tmp[dim]           = (tmp[dim] - 1 + dd->nc[dim]) % dd->nc[dim];
1292         dd->neighbor[d][1] = ddcoord2ddnodeid(dd, tmp);
1293         if (debug)
1294         {
1295             fprintf(debug, "DD rank %d neighbor ranks in dir %d are + %d - %d\n",
1296                     dd->rank, dim,
1297                     dd->neighbor[d][0],
1298                     dd->neighbor[d][1]);
1299         }
1300     }
1301
1302     int nzone  = (1 << dd->ndim);
1303     GMX_ASSERT(dd->ndim < DIM, "Invalid number of dimensions");
1304     int nizone = (1 << std::max(dd->ndim - 1, 0));
1305     assert(nizone >= 1 && nizone <= DD_MAXIZONE);
1306
1307     zones = &dd->comm->zones;
1308
1309     for (i = 0; i < nzone; i++)
1310     {
1311         m = 0;
1312         clear_ivec(zones->shift[i]);
1313         for (d = 0; d < dd->ndim; d++)
1314         {
1315             zones->shift[i][dd->dim[d]] = dd_zo[i][m++];
1316         }
1317     }
1318
1319     zones->n = nzone;
1320     for (i = 0; i < nzone; i++)
1321     {
1322         for (d = 0; d < DIM; d++)
1323         {
1324             s[d] = dd->ci[d] - zones->shift[i][d];
1325             if (s[d] < 0)
1326             {
1327                 s[d] += dd->nc[d];
1328             }
1329             else if (s[d] >= dd->nc[d])
1330             {
1331                 s[d] -= dd->nc[d];
1332             }
1333         }
1334     }
1335     zones->nizone = nizone;
1336     for (i = 0; i < zones->nizone; i++)
1337     {
1338         assert(ddNonbondedZonePairRanges[i][0] == i);
1339
1340         izone     = &zones->izone[i];
1341         /* dd_zp3 is for 3D decomposition, for fewer dimensions use only
1342          * j-zones up to nzone.
1343          */
1344         izone->j0 = std::min(ddNonbondedZonePairRanges[i][1], nzone);
1345         izone->j1 = std::min(ddNonbondedZonePairRanges[i][2], nzone);
1346         for (dim = 0; dim < DIM; dim++)
1347         {
1348             if (dd->nc[dim] == 1)
1349             {
1350                 /* All shifts should be allowed */
1351                 izone->shift0[dim] = -1;
1352                 izone->shift1[dim] = 1;
1353             }
1354             else
1355             {
1356                 /* Determine the min/max j-zone shift wrt the i-zone */
1357                 izone->shift0[dim] = 1;
1358                 izone->shift1[dim] = -1;
1359                 for (j = izone->j0; j < izone->j1; j++)
1360                 {
1361                     int shift_diff = zones->shift[j][dim] - zones->shift[i][dim];
1362                     if (shift_diff < izone->shift0[dim])
1363                     {
1364                         izone->shift0[dim] = shift_diff;
1365                     }
1366                     if (shift_diff > izone->shift1[dim])
1367                     {
1368                         izone->shift1[dim] = shift_diff;
1369                     }
1370                 }
1371             }
1372         }
1373     }
1374
1375     if (!isDlbDisabled(dd->comm))
1376     {
1377         dd->comm->cellsizesWithDlb.resize(dd->ndim);
1378     }
1379
1380     if (dd->comm->ddSettings.recordLoad)
1381     {
1382         make_load_communicators(dd);
1383     }
1384 }
1385
1386 static void make_pp_communicator(const gmx::MDLogger  &mdlog,
1387                                  gmx_domdec_t         *dd,
1388                                  t_commrec gmx_unused *cr,
1389                                  bool gmx_unused       reorder)
1390 {
1391 #if GMX_MPI
1392     gmx_domdec_comm_t *comm;
1393     int                rank, *buf;
1394     ivec               periods;
1395     MPI_Comm           comm_cart;
1396
1397     comm = dd->comm;
1398
1399     if (comm->bCartesianPP)
1400     {
1401         /* Set up cartesian communication for the particle-particle part */
1402         GMX_LOG(mdlog.info).appendTextFormatted(
1403                 "Will use a Cartesian communicator: %d x %d x %d",
1404                 dd->nc[XX], dd->nc[YY], dd->nc[ZZ]);
1405
1406         for (int i = 0; i < DIM; i++)
1407         {
1408             periods[i] = TRUE;
1409         }
1410         MPI_Cart_create(cr->mpi_comm_mygroup, DIM, dd->nc, periods, static_cast<int>(reorder),
1411                         &comm_cart);
1412         /* We overwrite the old communicator with the new cartesian one */
1413         cr->mpi_comm_mygroup = comm_cart;
1414     }
1415
1416     dd->mpi_comm_all = cr->mpi_comm_mygroup;
1417     MPI_Comm_rank(dd->mpi_comm_all, &dd->rank);
1418
1419     if (comm->bCartesianPP_PME)
1420     {
1421         /* Since we want to use the original cartesian setup for sim,
1422          * and not the one after split, we need to make an index.
1423          */
1424         snew(comm->ddindex2ddnodeid, dd->nnodes);
1425         comm->ddindex2ddnodeid[dd_index(dd->nc, dd->ci)] = dd->rank;
1426         gmx_sumi(dd->nnodes, comm->ddindex2ddnodeid, cr);
1427         /* Get the rank of the DD master,
1428          * above we made sure that the master node is a PP node.
1429          */
1430         if (MASTER(cr))
1431         {
1432             rank = dd->rank;
1433         }
1434         else
1435         {
1436             rank = 0;
1437         }
1438         MPI_Allreduce(&rank, &dd->masterrank, 1, MPI_INT, MPI_SUM, dd->mpi_comm_all);
1439     }
1440     else if (comm->bCartesianPP)
1441     {
1442         if (cr->npmenodes == 0)
1443         {
1444             /* The PP communicator is also
1445              * the communicator for this simulation
1446              */
1447             cr->mpi_comm_mysim = cr->mpi_comm_mygroup;
1448         }
1449         cr->nodeid = dd->rank;
1450
1451         MPI_Cart_coords(dd->mpi_comm_all, dd->rank, DIM, dd->ci);
1452
1453         /* We need to make an index to go from the coordinates
1454          * to the nodeid of this simulation.
1455          */
1456         snew(comm->ddindex2simnodeid, dd->nnodes);
1457         snew(buf, dd->nnodes);
1458         if (thisRankHasDuty(cr, DUTY_PP))
1459         {
1460             buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
1461         }
1462         /* Communicate the ddindex to simulation nodeid index */
1463         MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
1464                       cr->mpi_comm_mysim);
1465         sfree(buf);
1466
1467         /* Determine the master coordinates and rank.
1468          * The DD master should be the same node as the master of this sim.
1469          */
1470         for (int i = 0; i < dd->nnodes; i++)
1471         {
1472             if (comm->ddindex2simnodeid[i] == 0)
1473             {
1474                 ddindex2xyz(dd->nc, i, dd->master_ci);
1475                 MPI_Cart_rank(dd->mpi_comm_all, dd->master_ci, &dd->masterrank);
1476             }
1477         }
1478         if (debug)
1479         {
1480             fprintf(debug, "The master rank is %d\n", dd->masterrank);
1481         }
1482     }
1483     else
1484     {
1485         /* No Cartesian communicators */
1486         /* We use the rank in dd->comm->all as DD index */
1487         ddindex2xyz(dd->nc, dd->rank, dd->ci);
1488         /* The simulation master nodeid is 0, so the DD master rank is also 0 */
1489         dd->masterrank = 0;
1490         clear_ivec(dd->master_ci);
1491     }
1492 #endif
1493
1494     GMX_LOG(mdlog.info).appendTextFormatted(
1495             "Domain decomposition rank %d, coordinates %d %d %d\n",
1496             dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
1497     if (debug)
1498     {
1499         fprintf(debug,
1500                 "Domain decomposition rank %d, coordinates %d %d %d\n\n",
1501                 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
1502     }
1503 }
1504
1505 static void receive_ddindex2simnodeid(gmx_domdec_t         *dd,
1506                                       t_commrec            *cr)
1507 {
1508 #if GMX_MPI
1509     gmx_domdec_comm_t *comm = dd->comm;
1510
1511     if (!comm->bCartesianPP_PME && comm->bCartesianPP)
1512     {
1513         int *buf;
1514         snew(comm->ddindex2simnodeid, dd->nnodes);
1515         snew(buf, dd->nnodes);
1516         if (thisRankHasDuty(cr, DUTY_PP))
1517         {
1518             buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
1519         }
1520         /* Communicate the ddindex to simulation nodeid index */
1521         MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
1522                       cr->mpi_comm_mysim);
1523         sfree(buf);
1524     }
1525 #else
1526     GMX_UNUSED_VALUE(dd);
1527     GMX_UNUSED_VALUE(cr);
1528 #endif
1529 }
1530
1531 static void split_communicator(const gmx::MDLogger &mdlog,
1532                                t_commrec *cr, gmx_domdec_t *dd,
1533                                DdRankOrder gmx_unused rankOrder,
1534                                bool gmx_unused reorder)
1535 {
1536     gmx_domdec_comm_t *comm;
1537     int                i;
1538     gmx_bool           bDiv[DIM];
1539 #if GMX_MPI
1540     MPI_Comm           comm_cart;
1541 #endif
1542
1543     comm = dd->comm;
1544
1545     if (comm->bCartesianPP)
1546     {
1547         for (i = 1; i < DIM; i++)
1548         {
1549             bDiv[i] = ((cr->npmenodes*dd->nc[i]) % (dd->nnodes) == 0);
1550         }
1551         if (bDiv[YY] || bDiv[ZZ])
1552         {
1553             comm->bCartesianPP_PME = TRUE;
1554             /* If we have 2D PME decomposition, which is always in x+y,
1555              * we stack the PME only nodes in z.
1556              * Otherwise we choose the direction that provides the thinnest slab
1557              * of PME only nodes as this will have the least effect
1558              * on the PP communication.
1559              * But for the PME communication the opposite might be better.
1560              */
1561             if (bDiv[ZZ] && (comm->npmenodes_y > 1 ||
1562                              !bDiv[YY] ||
1563                              dd->nc[YY] > dd->nc[ZZ]))
1564             {
1565                 comm->cartpmedim = ZZ;
1566             }
1567             else
1568             {
1569                 comm->cartpmedim = YY;
1570             }
1571             comm->ntot[comm->cartpmedim]
1572                 += (cr->npmenodes*dd->nc[comm->cartpmedim])/dd->nnodes;
1573         }
1574         else
1575         {
1576             GMX_LOG(mdlog.info).appendTextFormatted(
1577                     "Number of PME-only ranks (%d) is not a multiple of nx*ny (%d*%d) or nx*nz (%d*%d)",
1578                     cr->npmenodes, dd->nc[XX], dd->nc[YY], dd->nc[XX], dd->nc[ZZ]);
1579             GMX_LOG(mdlog.info).appendText("Will not use a Cartesian communicator for PP <-> PME\n");
1580         }
1581     }
1582
1583     if (comm->bCartesianPP_PME)
1584     {
1585 #if GMX_MPI
1586         int  rank;
1587         ivec periods;
1588
1589         GMX_LOG(mdlog.info).appendTextFormatted(
1590                 "Will use a Cartesian communicator for PP <-> PME: %d x %d x %d",
1591                 comm->ntot[XX], comm->ntot[YY], comm->ntot[ZZ]);
1592
1593         for (i = 0; i < DIM; i++)
1594         {
1595             periods[i] = TRUE;
1596         }
1597         MPI_Cart_create(cr->mpi_comm_mysim, DIM, comm->ntot, periods, static_cast<int>(reorder),
1598                         &comm_cart);
1599         MPI_Comm_rank(comm_cart, &rank);
1600         if (MASTER(cr) && rank != 0)
1601         {
1602             gmx_fatal(FARGS, "MPI rank 0 was renumbered by MPI_Cart_create, we do not allow this");
1603         }
1604
1605         /* With this assigment we loose the link to the original communicator
1606          * which will usually be MPI_COMM_WORLD, unless have multisim.
1607          */
1608         cr->mpi_comm_mysim = comm_cart;
1609         cr->sim_nodeid     = rank;
1610
1611         MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, dd->ci);
1612
1613         GMX_LOG(mdlog.info).appendTextFormatted(
1614                 "Cartesian rank %d, coordinates %d %d %d\n",
1615                 cr->sim_nodeid, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
1616
1617         if (dd->ci[comm->cartpmedim] < dd->nc[comm->cartpmedim])
1618         {
1619             cr->duty = DUTY_PP;
1620         }
1621         if (cr->npmenodes == 0 ||
1622             dd->ci[comm->cartpmedim] >= dd->nc[comm->cartpmedim])
1623         {
1624             cr->duty = DUTY_PME;
1625         }
1626
1627         /* Split the sim communicator into PP and PME only nodes */
1628         MPI_Comm_split(cr->mpi_comm_mysim,
1629                        getThisRankDuties(cr),
1630                        dd_index(comm->ntot, dd->ci),
1631                        &cr->mpi_comm_mygroup);
1632 #endif
1633     }
1634     else
1635     {
1636         switch (rankOrder)
1637         {
1638             case DdRankOrder::pp_pme:
1639                 GMX_LOG(mdlog.info).appendText("Order of the ranks: PP first, PME last");
1640                 break;
1641             case DdRankOrder::interleave:
1642                 /* Interleave the PP-only and PME-only ranks */
1643                 GMX_LOG(mdlog.info).appendText("Interleaving PP and PME ranks");
1644                 comm->pmenodes = dd_interleaved_pme_ranks(dd);
1645                 break;
1646             case DdRankOrder::cartesian:
1647                 break;
1648             default:
1649                 gmx_fatal(FARGS, "Invalid ddRankOrder=%d", static_cast<int>(rankOrder));
1650         }
1651
1652         if (dd_simnode2pmenode(dd, cr, cr->sim_nodeid) == -1)
1653         {
1654             cr->duty = DUTY_PME;
1655         }
1656         else
1657         {
1658             cr->duty = DUTY_PP;
1659         }
1660 #if GMX_MPI
1661         /* Split the sim communicator into PP and PME only nodes */
1662         MPI_Comm_split(cr->mpi_comm_mysim,
1663                        getThisRankDuties(cr),
1664                        cr->nodeid,
1665                        &cr->mpi_comm_mygroup);
1666         MPI_Comm_rank(cr->mpi_comm_mygroup, &cr->nodeid);
1667 #endif
1668     }
1669
1670     GMX_LOG(mdlog.info).appendTextFormatted(
1671             "This rank does only %s work.\n",
1672             thisRankHasDuty(cr, DUTY_PP) ? "particle-particle" : "PME-mesh");
1673 }
1674
1675 /*! \brief Generates the MPI communicators for domain decomposition */
1676 static void make_dd_communicators(const gmx::MDLogger &mdlog,
1677                                   t_commrec *cr,
1678                                   gmx_domdec_t *dd, DdRankOrder ddRankOrder)
1679 {
1680     gmx_domdec_comm_t *comm;
1681     bool               CartReorder;
1682
1683     comm = dd->comm;
1684
1685     copy_ivec(dd->nc, comm->ntot);
1686
1687     comm->bCartesianPP     = (ddRankOrder == DdRankOrder::cartesian);
1688     comm->bCartesianPP_PME = FALSE;
1689
1690     /* Reorder the nodes by default. This might change the MPI ranks.
1691      * Real reordering is only supported on very few architectures,
1692      * Blue Gene is one of them.
1693      */
1694     CartReorder = getenv("GMX_NO_CART_REORDER") == nullptr;
1695
1696     if (cr->npmenodes > 0)
1697     {
1698         /* Split the communicator into a PP and PME part */
1699         split_communicator(mdlog, cr, dd, ddRankOrder, CartReorder);
1700         if (comm->bCartesianPP_PME)
1701         {
1702             /* We (possibly) reordered the nodes in split_communicator,
1703              * so it is no longer required in make_pp_communicator.
1704              */
1705             CartReorder = FALSE;
1706         }
1707     }
1708     else
1709     {
1710         /* All nodes do PP and PME */
1711 #if GMX_MPI
1712         /* We do not require separate communicators */
1713         cr->mpi_comm_mygroup = cr->mpi_comm_mysim;
1714 #endif
1715     }
1716
1717     if (thisRankHasDuty(cr, DUTY_PP))
1718     {
1719         /* Copy or make a new PP communicator */
1720         make_pp_communicator(mdlog, dd, cr, CartReorder);
1721     }
1722     else
1723     {
1724         receive_ddindex2simnodeid(dd, cr);
1725     }
1726
1727     if (!thisRankHasDuty(cr, DUTY_PME))
1728     {
1729         /* Set up the commnuication to our PME node */
1730         dd->pme_nodeid           = dd_simnode2pmenode(dd, cr, cr->sim_nodeid);
1731         dd->pme_receive_vir_ener = receive_vir_ener(dd, cr);
1732         if (debug)
1733         {
1734             fprintf(debug, "My pme_nodeid %d receive ener %s\n",
1735                     dd->pme_nodeid, gmx::boolToString(dd->pme_receive_vir_ener));
1736         }
1737     }
1738     else
1739     {
1740         dd->pme_nodeid = -1;
1741     }
1742
1743     /* We can not use DDMASTER(dd), because dd->masterrank is set later */
1744     if (MASTER(cr))
1745     {
1746         dd->ma = std::make_unique<AtomDistribution>(dd->nc,
1747                                                     comm->cgs_gl.nr,
1748                                                     comm->cgs_gl.index[comm->cgs_gl.nr]);
1749     }
1750 }
1751
1752 static real *get_slb_frac(const gmx::MDLogger &mdlog,
1753                           const char *dir, int nc, const char *size_string)
1754 {
1755     real  *slb_frac, tot;
1756     int    i, n;
1757     double dbl;
1758
1759     slb_frac = nullptr;
1760     if (nc > 1 && size_string != nullptr)
1761     {
1762         GMX_LOG(mdlog.info).appendTextFormatted(
1763                 "Using static load balancing for the %s direction", dir);
1764         snew(slb_frac, nc);
1765         tot = 0;
1766         for (i = 0; i < nc; i++)
1767         {
1768             dbl = 0;
1769             sscanf(size_string, "%20lf%n", &dbl, &n);
1770             if (dbl == 0)
1771             {
1772                 gmx_fatal(FARGS, "Incorrect or not enough DD cell size entries for direction %s: '%s'", dir, size_string);
1773             }
1774             slb_frac[i]  = dbl;
1775             size_string += n;
1776             tot         += slb_frac[i];
1777         }
1778         /* Normalize */
1779         std::string relativeCellSizes = "Relative cell sizes:";
1780         for (i = 0; i < nc; i++)
1781         {
1782             slb_frac[i]       /= tot;
1783             relativeCellSizes += gmx::formatString(" %5.3f", slb_frac[i]);
1784         }
1785         GMX_LOG(mdlog.info).appendText(relativeCellSizes);
1786     }
1787
1788     return slb_frac;
1789 }
1790
1791 static int multi_body_bondeds_count(const gmx_mtop_t *mtop)
1792 {
1793     int                  n     = 0;
1794     gmx_mtop_ilistloop_t iloop = gmx_mtop_ilistloop_init(mtop);
1795     int                  nmol;
1796     while (const InteractionLists *ilists = gmx_mtop_ilistloop_next(iloop, &nmol))
1797     {
1798         for (auto &ilist : extractILists(*ilists, IF_BOND))
1799         {
1800             if (NRAL(ilist.functionType) >  2)
1801             {
1802                 n += nmol*(ilist.iatoms.size()/ilistStride(ilist));
1803             }
1804         }
1805     }
1806
1807     return n;
1808 }
1809
1810 static int dd_getenv(const gmx::MDLogger &mdlog,
1811                      const char *env_var, int def)
1812 {
1813     char *val;
1814     int   nst;
1815
1816     nst = def;
1817     val = getenv(env_var);
1818     if (val)
1819     {
1820         if (sscanf(val, "%20d", &nst) <= 0)
1821         {
1822             nst = 1;
1823         }
1824         GMX_LOG(mdlog.info).appendTextFormatted(
1825                 "Found env.var. %s = %s, using value %d",
1826                 env_var, val, nst);
1827     }
1828
1829     return nst;
1830 }
1831
1832 static void check_dd_restrictions(const gmx_domdec_t  *dd,
1833                                   const t_inputrec    *ir,
1834                                   const gmx::MDLogger &mdlog)
1835 {
1836     if (ir->ePBC == epbcSCREW &&
1837         (dd->nc[XX] == 1 || dd->nc[YY] > 1 || dd->nc[ZZ] > 1))
1838     {
1839         gmx_fatal(FARGS, "With pbc=%s can only do domain decomposition in the x-direction", epbc_names[ir->ePBC]);
1840     }
1841
1842     if (ir->ns_type == ensSIMPLE)
1843     {
1844         gmx_fatal(FARGS, "Domain decomposition does not support simple neighbor searching, use grid searching or run with one MPI rank");
1845     }
1846
1847     if (ir->nstlist == 0)
1848     {
1849         gmx_fatal(FARGS, "Domain decomposition does not work with nstlist=0");
1850     }
1851
1852     if (ir->comm_mode == ecmANGULAR && ir->ePBC != epbcNONE)
1853     {
1854         GMX_LOG(mdlog.warning).appendText("comm-mode angular will give incorrect results when the comm group partially crosses a periodic boundary");
1855     }
1856 }
1857
1858 static real average_cellsize_min(const gmx_ddbox_t &ddbox,
1859                                  const ivec         numDomains)
1860 {
1861     real r = ddbox.box_size[XX];
1862     for (int d = 0; d < DIM; d++)
1863     {
1864         if (numDomains[d] > 1)
1865         {
1866             /* Check using the initial average cell size */
1867             r = std::min(r, ddbox.box_size[d]*ddbox.skew_fac[d]/numDomains[d]);
1868         }
1869     }
1870
1871     return r;
1872 }
1873
1874 /*! \brief Depending on the DLB initial value return the DLB switched off state or issue an error.
1875  */
1876 static DlbState forceDlbOffOrBail(DlbState             cmdlineDlbState,
1877                                   const std::string   &reasonStr,
1878                                   const gmx::MDLogger &mdlog)
1879 {
1880     std::string dlbNotSupportedErr  = "Dynamic load balancing requested, but ";
1881     std::string dlbDisableNote      = "NOTE: disabling dynamic load balancing as ";
1882
1883     if (cmdlineDlbState == DlbState::onUser)
1884     {
1885         gmx_fatal(FARGS, "%s", (dlbNotSupportedErr + reasonStr).c_str());
1886     }
1887     else if (cmdlineDlbState == DlbState::offCanTurnOn)
1888     {
1889         GMX_LOG(mdlog.info).appendText(dlbDisableNote + reasonStr);
1890     }
1891     return DlbState::offForever;
1892 }
1893
1894 /*! \brief Return the dynamic load balancer's initial state based on initial conditions and user inputs.
1895  *
1896  * This function parses the parameters of "-dlb" command line option setting
1897  * corresponding state values. Then it checks the consistency of the determined
1898  * state with other run parameters and settings. As a result, the initial state
1899  * may be altered or an error may be thrown if incompatibility of options is detected.
1900  *
1901  * \param [in] mdlog       Logger.
1902  * \param [in] dlbOption   Enum value for the DLB option.
1903  * \param [in] bRecordLoad True if the load balancer is recording load information.
1904  * \param [in] mdrunOptions  Options for mdrun.
1905  * \param [in] ir          Pointer mdrun to input parameters.
1906  * \returns                DLB initial/startup state.
1907  */
1908 static DlbState determineInitialDlbState(const gmx::MDLogger &mdlog,
1909                                          DlbOption dlbOption, gmx_bool bRecordLoad,
1910                                          const gmx::MdrunOptions &mdrunOptions,
1911                                          const t_inputrec *ir)
1912 {
1913     DlbState dlbState = DlbState::offCanTurnOn;
1914
1915     switch (dlbOption)
1916     {
1917         case DlbOption::turnOnWhenUseful: dlbState = DlbState::offCanTurnOn; break;
1918         case DlbOption::no:               dlbState = DlbState::offUser;      break;
1919         case DlbOption::yes:              dlbState = DlbState::onUser;       break;
1920         default: gmx_incons("Invalid dlbOption enum value");
1921     }
1922
1923     /* Reruns don't support DLB: bail or override auto mode */
1924     if (mdrunOptions.rerun)
1925     {
1926         std::string reasonStr = "it is not supported in reruns.";
1927         return forceDlbOffOrBail(dlbState, reasonStr, mdlog);
1928     }
1929
1930     /* Unsupported integrators */
1931     if (!EI_DYNAMICS(ir->eI))
1932     {
1933         auto reasonStr = gmx::formatString("it is only supported with dynamics, not with integrator '%s'.", EI(ir->eI));
1934         return forceDlbOffOrBail(dlbState, reasonStr, mdlog);
1935     }
1936
1937     /* Without cycle counters we can't time work to balance on */
1938     if (!bRecordLoad)
1939     {
1940         std::string reasonStr = "cycle counters unsupported or not enabled in the operating system kernel.";
1941         return forceDlbOffOrBail(dlbState, reasonStr, mdlog);
1942     }
1943
1944     if (mdrunOptions.reproducible)
1945     {
1946         std::string reasonStr = "you started a reproducible run.";
1947         switch (dlbState)
1948         {
1949             case DlbState::offUser:
1950                 break;
1951             case DlbState::offForever:
1952                 GMX_RELEASE_ASSERT(false, "DlbState::offForever is not a valid initial state");
1953                 break;
1954             case DlbState::offCanTurnOn:
1955                 return forceDlbOffOrBail(dlbState, reasonStr, mdlog);
1956             case DlbState::onCanTurnOff:
1957                 GMX_RELEASE_ASSERT(false, "DlbState::offCanTurnOff is not a valid initial state");
1958                 break;
1959             case DlbState::onUser:
1960                 return forceDlbOffOrBail(dlbState, reasonStr + " In load balanced runs binary reproducibility cannot be ensured.", mdlog);
1961             default:
1962                 gmx_fatal(FARGS, "Death horror: undefined case (%d) for load balancing choice", static_cast<int>(dlbState));
1963         }
1964     }
1965
1966     return dlbState;
1967 }
1968
1969 static void set_dd_dim(const gmx::MDLogger &mdlog, gmx_domdec_t *dd)
1970 {
1971     dd->ndim = 0;
1972     if (getenv("GMX_DD_ORDER_ZYX") != nullptr)
1973     {
1974         /* Decomposition order z,y,x */
1975         GMX_LOG(mdlog.info).appendText("Using domain decomposition order z, y, x");
1976         for (int dim = DIM-1; dim >= 0; dim--)
1977         {
1978             if (dd->nc[dim] > 1)
1979             {
1980                 dd->dim[dd->ndim++] = dim;
1981             }
1982         }
1983     }
1984     else
1985     {
1986         /* Decomposition order x,y,z */
1987         for (int dim = 0; dim < DIM; dim++)
1988         {
1989             if (dd->nc[dim] > 1)
1990             {
1991                 dd->dim[dd->ndim++] = dim;
1992             }
1993         }
1994     }
1995
1996     if (dd->ndim == 0)
1997     {
1998         /* Set dim[0] to avoid extra checks on ndim in several places */
1999         dd->dim[0] = XX;
2000     }
2001 }
2002
2003 static gmx_domdec_comm_t *init_dd_comm()
2004 {
2005     gmx_domdec_comm_t *comm = new gmx_domdec_comm_t;
2006
2007     comm->n_load_have      = 0;
2008     comm->n_load_collect   = 0;
2009
2010     comm->haveTurnedOffDlb = false;
2011
2012     for (int i = 0; i < static_cast<int>(DDAtomRanges::Type::Number); i++)
2013     {
2014         comm->sum_nat[i] = 0;
2015     }
2016     comm->ndecomp   = 0;
2017     comm->nload     = 0;
2018     comm->load_step = 0;
2019     comm->load_sum  = 0;
2020     comm->load_max  = 0;
2021     clear_ivec(comm->load_lim);
2022     comm->load_mdf  = 0;
2023     comm->load_pme  = 0;
2024
2025     /* This should be replaced by a unique pointer */
2026     comm->balanceRegion = ddBalanceRegionAllocate();
2027
2028     return comm;
2029 }
2030
2031 /* Returns whether mtop contains constraints and/or vsites */
2032 static bool systemHasConstraintsOrVsites(const gmx_mtop_t &mtop)
2033 {
2034     auto ilistLoop = gmx_mtop_ilistloop_init(mtop);
2035     int  nmol;
2036     while (const InteractionLists *ilists = gmx_mtop_ilistloop_next(ilistLoop, &nmol))
2037     {
2038         if (!extractILists(*ilists, IF_CONSTRAINT | IF_VSITE).empty())
2039         {
2040             return true;
2041         }
2042     }
2043
2044     return false;
2045 }
2046
2047 static void setupUpdateGroups(const gmx::MDLogger &mdlog,
2048                               const gmx_mtop_t    &mtop,
2049                               const t_inputrec    &inputrec,
2050                               const real           cutoffMargin,
2051                               DDSystemInfo        *systemInfo)
2052 {
2053     /* When we have constraints and/or vsites, it is beneficial to use
2054      * update groups (when possible) to allow independent update of groups.
2055      */
2056     if (!systemHasConstraintsOrVsites(mtop))
2057     {
2058         /* No constraints or vsites, atoms can be updated independently */
2059         return;
2060     }
2061
2062     systemInfo->updateGroupingPerMoleculetype = gmx::makeUpdateGroups(mtop);
2063     systemInfo->useUpdateGroups               =
2064         (!systemInfo->updateGroupingPerMoleculetype.empty() &&
2065          getenv("GMX_NO_UPDATEGROUPS") == nullptr);
2066
2067     if (systemInfo->useUpdateGroups)
2068     {
2069         int numUpdateGroups = 0;
2070         for (const auto &molblock : mtop.molblock)
2071         {
2072             numUpdateGroups += molblock.nmol*systemInfo->updateGroupingPerMoleculetype[molblock.type].numBlocks();
2073         }
2074
2075         systemInfo->maxUpdateGroupRadius =
2076             computeMaxUpdateGroupRadius(mtop,
2077                                         systemInfo->updateGroupingPerMoleculetype,
2078                                         maxReferenceTemperature(inputrec));
2079
2080         /* To use update groups, the large domain-to-domain cutoff distance
2081          * should be compatible with the box size.
2082          */
2083         systemInfo->useUpdateGroups = (atomToAtomIntoDomainToDomainCutoff(*systemInfo, 0) < cutoffMargin);
2084
2085         if (systemInfo->useUpdateGroups)
2086         {
2087             GMX_LOG(mdlog.info).appendTextFormatted(
2088                     "Using update groups, nr %d, average size %.1f atoms, max. radius %.3f nm\n",
2089                     numUpdateGroups,
2090                     mtop.natoms/static_cast<double>(numUpdateGroups),
2091                     systemInfo->maxUpdateGroupRadius);
2092         }
2093         else
2094         {
2095             GMX_LOG(mdlog.info).appendTextFormatted("The combination of rlist and box size prohibits the use of update groups\n");
2096             systemInfo->updateGroupingPerMoleculetype.clear();
2097         }
2098     }
2099 }
2100
2101 UnitCellInfo::UnitCellInfo(const t_inputrec &ir) :
2102     npbcdim(ePBC2npbcdim(ir.ePBC)),
2103     numBoundedDimensions(inputrec2nboundeddim(&ir)),
2104     ddBoxIsDynamic(numBoundedDimensions < DIM || inputrecDynamicBox(&ir)),
2105     haveScrewPBC(ir.ePBC == epbcSCREW)
2106 {
2107 }
2108
2109 /*! \brief Generate the simulation system information */
2110 static DDSystemInfo
2111 getSystemInfo(const gmx::MDLogger           &mdlog,
2112               t_commrec                     *cr,
2113               const DomdecOptions           &options,
2114               const gmx_mtop_t              *mtop,
2115               const t_inputrec              *ir,
2116               const matrix                   box,
2117               gmx::ArrayRef<const gmx::RVec> xGlobal)
2118 {
2119     const real         tenPercentMargin = 1.1;
2120
2121     DDSystemInfo       systemInfo;
2122
2123     /* We need to decide on update groups early, as this affects communication distances */
2124     systemInfo.useUpdateGroups = false;
2125     if (ir->cutoff_scheme == ecutsVERLET)
2126     {
2127         real cutoffMargin = std::sqrt(max_cutoff2(ir->ePBC, box)) - ir->rlist;
2128         setupUpdateGroups(mdlog, *mtop, *ir, cutoffMargin, &systemInfo);
2129     }
2130
2131     // TODO: Check whether all bondeds are within update groups
2132     systemInfo.haveInterDomainBondeds          = (mtop->natoms > gmx_mtop_num_molecules(*mtop) ||
2133                                                   mtop->bIntermolecularInteractions);
2134     systemInfo.haveInterDomainMultiBodyBondeds = (multi_body_bondeds_count(mtop) > 0);
2135
2136     if (systemInfo.useUpdateGroups)
2137     {
2138         systemInfo.haveSplitConstraints = false;
2139         systemInfo.haveSplitSettles     = false;
2140     }
2141     else
2142     {
2143         systemInfo.haveSplitConstraints = gmx::inter_charge_group_constraints(*mtop);
2144         systemInfo.haveSplitSettles     = gmx::inter_charge_group_settles(*mtop);
2145     }
2146
2147     if (ir->rlist == 0)
2148     {
2149         /* Set the cut-off to some very large value,
2150          * so we don't need if statements everywhere in the code.
2151          * We use sqrt, since the cut-off is squared in some places.
2152          */
2153         systemInfo.cutoff = GMX_CUTOFF_INF;
2154     }
2155     else
2156     {
2157         systemInfo.cutoff = atomToAtomIntoDomainToDomainCutoff(systemInfo, ir->rlist);
2158     }
2159     systemInfo.minCutoffForMultiBody = 0;
2160
2161     /* Determine the minimum cell size limit, affected by many factors */
2162     systemInfo.cellsizeLimit             = 0;
2163     systemInfo.filterBondedCommunication = false;
2164
2165     /* We do not allow home atoms to move beyond the neighboring domain
2166      * between domain decomposition steps, which limits the cell size.
2167      * Get an estimate of cell size limit due to atom displacement.
2168      * In most cases this is a large overestimate, because it assumes
2169      * non-interaction atoms.
2170      * We set the chance to 1 in a trillion steps.
2171      */
2172     constexpr real c_chanceThatAtomMovesBeyondDomain = 1e-12;
2173     const real     limitForAtomDisplacement          =
2174         minCellSizeForAtomDisplacement(*mtop, *ir,
2175                                        systemInfo.updateGroupingPerMoleculetype,
2176                                        c_chanceThatAtomMovesBeyondDomain);
2177     GMX_LOG(mdlog.info).appendTextFormatted(
2178             "Minimum cell size due to atom displacement: %.3f nm",
2179             limitForAtomDisplacement);
2180
2181     systemInfo.cellsizeLimit = std::max(systemInfo.cellsizeLimit,
2182                                         limitForAtomDisplacement);
2183
2184     /* TODO: PME decomposition currently requires atoms not to be more than
2185      *       2/3 of comm->cutoff, which is >=rlist, outside of their domain.
2186      *       In nearly all cases, limitForAtomDisplacement will be smaller
2187      *       than 2/3*rlist, so the PME requirement is satisfied.
2188      *       But it would be better for both correctness and performance
2189      *       to use limitForAtomDisplacement instead of 2/3*comm->cutoff.
2190      *       Note that we would need to improve the pairlist buffer case.
2191      */
2192
2193     if (systemInfo.haveInterDomainBondeds)
2194     {
2195         if (options.minimumCommunicationRange > 0)
2196         {
2197             systemInfo.minCutoffForMultiBody =
2198                 atomToAtomIntoDomainToDomainCutoff(systemInfo, options.minimumCommunicationRange);
2199             if (options.useBondedCommunication)
2200             {
2201                 systemInfo.filterBondedCommunication = (systemInfo.minCutoffForMultiBody > systemInfo.cutoff);
2202             }
2203             else
2204             {
2205                 systemInfo.cutoff = std::max(systemInfo.cutoff,
2206                                              systemInfo.minCutoffForMultiBody);
2207             }
2208         }
2209         else if (ir->bPeriodicMols)
2210         {
2211             /* Can not easily determine the required cut-off */
2212             GMX_LOG(mdlog.warning).appendText("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.");
2213             systemInfo.minCutoffForMultiBody = systemInfo.cutoff/2;
2214         }
2215         else
2216         {
2217             real r_2b, r_mb;
2218
2219             if (MASTER(cr))
2220             {
2221                 dd_bonded_cg_distance(mdlog, mtop, ir, as_rvec_array(xGlobal.data()), box,
2222                                       options.checkBondedInteractions,
2223                                       &r_2b, &r_mb);
2224             }
2225             gmx_bcast(sizeof(r_2b), &r_2b, cr);
2226             gmx_bcast(sizeof(r_mb), &r_mb, cr);
2227
2228             /* We use an initial margin of 10% for the minimum cell size,
2229              * except when we are just below the non-bonded cut-off.
2230              */
2231             if (options.useBondedCommunication)
2232             {
2233                 if (std::max(r_2b, r_mb) > systemInfo.cutoff)
2234                 {
2235                     const real r_bonded              = std::max(r_2b, r_mb);
2236                     systemInfo.minCutoffForMultiBody = tenPercentMargin*r_bonded;
2237                     /* This is the (only) place where we turn on the filtering */
2238                     systemInfo.filterBondedCommunication = true;
2239                 }
2240                 else
2241                 {
2242                     const real r_bonded              = r_mb;
2243                     systemInfo.minCutoffForMultiBody = std::min(tenPercentMargin*r_bonded,
2244                                                                 systemInfo.cutoff);
2245                 }
2246                 /* We determine cutoff_mbody later */
2247                 systemInfo.increaseMultiBodyCutoff = true;
2248             }
2249             else
2250             {
2251                 /* No special bonded communication,
2252                  * simply increase the DD cut-off.
2253                  */
2254                 systemInfo.minCutoffForMultiBody = tenPercentMargin*std::max(r_2b, r_mb);
2255                 systemInfo.cutoff                = std::max(systemInfo.cutoff,
2256                                                             systemInfo.minCutoffForMultiBody);
2257             }
2258         }
2259         GMX_LOG(mdlog.info).appendTextFormatted(
2260                 "Minimum cell size due to bonded interactions: %.3f nm",
2261                 systemInfo.minCutoffForMultiBody);
2262
2263         systemInfo.cellsizeLimit = std::max(systemInfo.cellsizeLimit,
2264                                             systemInfo.minCutoffForMultiBody);
2265     }
2266
2267     systemInfo.constraintCommunicationRange = 0;
2268     if (systemInfo.haveSplitConstraints && options.constraintCommunicationRange <= 0)
2269     {
2270         /* There is a cell size limit due to the constraints (P-LINCS) */
2271         systemInfo.constraintCommunicationRange = gmx::constr_r_max(mdlog, mtop, ir);
2272         GMX_LOG(mdlog.info).appendTextFormatted(
2273                 "Estimated maximum distance required for P-LINCS: %.3f nm",
2274                 systemInfo.constraintCommunicationRange);
2275         if (systemInfo.constraintCommunicationRange > systemInfo.cellsizeLimit)
2276         {
2277             GMX_LOG(mdlog.info).appendText("This distance will limit the DD cell size, you can override this with -rcon");
2278         }
2279     }
2280     else if (options.constraintCommunicationRange > 0)
2281     {
2282         /* Here we do not check for dd->splitConstraints.
2283          * because one can also set a cell size limit for virtual sites only
2284          * and at this point we don't know yet if there are intercg v-sites.
2285          */
2286         GMX_LOG(mdlog.info).appendTextFormatted(
2287                 "User supplied maximum distance required for P-LINCS: %.3f nm",
2288                 options.constraintCommunicationRange);
2289         systemInfo.constraintCommunicationRange = options.constraintCommunicationRange;
2290     }
2291     systemInfo.cellsizeLimit = std::max(systemInfo.cellsizeLimit,
2292                                         systemInfo.constraintCommunicationRange);
2293
2294     return systemInfo;
2295 }
2296
2297 /*! \brief Set the cell size and interaction limits, as well as the DD grid */
2298 static void set_dd_limits_and_grid(const gmx::MDLogger &mdlog,
2299                                    t_commrec *cr, gmx_domdec_t *dd,
2300                                    const DomdecOptions &options,
2301                                    const DDSettings &ddSettings,
2302                                    const gmx_mtop_t *mtop,
2303                                    const t_inputrec *ir,
2304                                    const matrix box,
2305                                    gmx::ArrayRef<const gmx::RVec> xGlobal,
2306                                    gmx_ddbox_t *ddbox)
2307 {
2308     gmx_domdec_comm_t *comm = dd->comm;
2309     comm->ddSettings        = ddSettings;
2310
2311     /* Initialize to GPU share count to 0, might change later */
2312     comm->nrank_gpu_shared = 0;
2313
2314     comm->dlbState         = comm->ddSettings.initialDlbState;
2315     dd_dlb_set_should_check_whether_to_turn_dlb_on(dd, TRUE);
2316     /* To consider turning DLB on after 2*nstlist steps we need to check
2317      * at partitioning count 3. Thus we need to increase the first count by 2.
2318      */
2319     comm->ddPartioningCountFirstDlbOff += 2;
2320
2321     GMX_LOG(mdlog.info).appendTextFormatted(
2322             "Dynamic load balancing: %s", edlbs_names[int(comm->dlbState)]);
2323
2324     comm->bPMELoadBalDLBLimits = FALSE;
2325
2326     /* Allocate the charge group/atom sorting struct */
2327     comm->sort = std::make_unique<gmx_domdec_sort_t>();
2328
2329     /* Generate the simulation system information */
2330     comm->systemInfo               = getSystemInfo(mdlog, cr, options, mtop, ir, box, xGlobal);
2331     const DDSystemInfo &systemInfo = comm->systemInfo;
2332
2333     if (systemInfo.useUpdateGroups)
2334     {
2335         /* Note: We would like to use dd->nnodes for the atom count estimate,
2336          *       but that is not yet available here. But this anyhow only
2337          *       affect performance up to the second dd_partition_system call.
2338          */
2339         const int homeAtomCountEstimate =  mtop->natoms/cr->nnodes;
2340         comm->updateGroupsCog =
2341             std::make_unique<gmx::UpdateGroupsCog>(*mtop,
2342                                                    systemInfo.updateGroupingPerMoleculetype,
2343                                                    maxReferenceTemperature(*ir),
2344                                                    homeAtomCountEstimate);
2345     }
2346
2347     comm->cgs_gl = gmx_mtop_global_cgs(mtop);
2348
2349     DDSetup ddSetup;
2350     if (options.numCells[XX] > 0)
2351     {
2352         copy_ivec(options.numCells, ddSetup.numDomains);
2353         set_ddbox_cr(*cr, &ddSetup.numDomains, *ir, box, xGlobal, ddbox);
2354
2355         if (options.numPmeRanks >= 0)
2356         {
2357             ddSetup.numPmeRanks = options.numPmeRanks;
2358         }
2359         else
2360         {
2361             /* When the DD grid is set explicitly and -npme is set to auto,
2362              * don't use PME ranks. We check later if the DD grid is
2363              * compatible with the total number of ranks.
2364              */
2365             ddSetup.numPmeRanks = 0;
2366         }
2367     }
2368     else
2369     {
2370         set_ddbox_cr(*cr, nullptr, *ir, box, xGlobal, ddbox);
2371
2372         /* We need to choose the optimal DD grid and possibly PME nodes */
2373         ddSetup =
2374             dd_choose_grid(mdlog, cr, ir, mtop, box, ddbox,
2375                            options.numPmeRanks,
2376                            !isDlbDisabled(comm),
2377                            options.dlbScaling,
2378                            systemInfo);
2379
2380         if (ddSetup.numDomains[XX] == 0)
2381         {
2382             char     buf[STRLEN];
2383             gmx_bool bC = (systemInfo.haveSplitConstraints &&
2384                            systemInfo.constraintCommunicationRange > systemInfo.minCutoffForMultiBody);
2385             sprintf(buf, "Change the number of ranks or mdrun option %s%s%s",
2386                     !bC ? "-rdd" : "-rcon",
2387                     comm->dlbState != DlbState::offUser ? " or -dds" : "",
2388                     bC ? " or your LINCS settings" : "");
2389
2390             gmx_fatal_collective(FARGS, cr->mpi_comm_mysim, MASTER(cr),
2391                                  "There is no domain decomposition for %d ranks that is compatible with the given box and a minimum cell size of %g nm\n"
2392                                  "%s\n"
2393                                  "Look in the log file for details on the domain decomposition",
2394                                  cr->nnodes - ddSetup.numPmeRanks,
2395                                  ddSetup.cellsizeLimit,
2396                                  buf);
2397         }
2398     }
2399
2400     const real acs = average_cellsize_min(*ddbox, ddSetup.numDomains);
2401     if (acs < systemInfo.cellsizeLimit)
2402     {
2403         if (options.numCells[XX] <= 0)
2404         {
2405             GMX_RELEASE_ASSERT(false, "dd_choose_grid() should return a grid that satisfies the cell size limits");
2406         }
2407         else
2408         {
2409             gmx_fatal_collective(FARGS, cr->mpi_comm_mysim, MASTER(cr),
2410                                  "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",
2411                                  acs, systemInfo.cellsizeLimit);
2412         }
2413     }
2414
2415     /* Set the DD setup given by ddSetup */
2416     cr->npmenodes = ddSetup.numPmeRanks;
2417     copy_ivec(ddSetup.numDomains, dd->nc);
2418     set_dd_dim(mdlog, dd);
2419
2420     GMX_LOG(mdlog.info).appendTextFormatted(
2421             "Domain decomposition grid %d x %d x %d, separate PME ranks %d",
2422             dd->nc[XX], dd->nc[YY], dd->nc[ZZ], cr->npmenodes);
2423
2424     dd->nnodes = dd->nc[XX]*dd->nc[YY]*dd->nc[ZZ];
2425     if (cr->nnodes - dd->nnodes != cr->npmenodes)
2426     {
2427         gmx_fatal_collective(FARGS, cr->mpi_comm_mysim, MASTER(cr),
2428                              "The size of the domain decomposition grid (%d) does not match the number of ranks (%d). The total number of ranks is %d",
2429                              dd->nnodes, cr->nnodes - cr->npmenodes, cr->nnodes);
2430     }
2431     if (cr->npmenodes > dd->nnodes)
2432     {
2433         gmx_fatal_collective(FARGS, cr->mpi_comm_mysim, MASTER(cr),
2434                              "The number of separate PME ranks (%d) is larger than the number of PP ranks (%d), this is not supported.", cr->npmenodes, dd->nnodes);
2435     }
2436     if (cr->npmenodes > 0)
2437     {
2438         comm->npmenodes = cr->npmenodes;
2439     }
2440     else
2441     {
2442         comm->npmenodes = dd->nnodes;
2443     }
2444
2445     if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
2446     {
2447         /* The following choices should match those
2448          * in comm_cost_est in domdec_setup.c.
2449          * Note that here the checks have to take into account
2450          * that the decomposition might occur in a different order than xyz
2451          * (for instance through the env.var. GMX_DD_ORDER_ZYX),
2452          * in which case they will not match those in comm_cost_est,
2453          * but since that is mainly for testing purposes that's fine.
2454          */
2455         if (dd->ndim >= 2 && dd->dim[0] == XX && dd->dim[1] == YY &&
2456             comm->npmenodes > dd->nc[XX] && comm->npmenodes % dd->nc[XX] == 0 &&
2457             getenv("GMX_PMEONEDD") == nullptr)
2458         {
2459             comm->npmedecompdim = 2;
2460             comm->npmenodes_x   = dd->nc[XX];
2461             comm->npmenodes_y   = comm->npmenodes/comm->npmenodes_x;
2462         }
2463         else
2464         {
2465             /* In case nc is 1 in both x and y we could still choose to
2466              * decompose pme in y instead of x, but we use x for simplicity.
2467              */
2468             comm->npmedecompdim = 1;
2469             if (dd->dim[0] == YY)
2470             {
2471                 comm->npmenodes_x = 1;
2472                 comm->npmenodes_y = comm->npmenodes;
2473             }
2474             else
2475             {
2476                 comm->npmenodes_x = comm->npmenodes;
2477                 comm->npmenodes_y = 1;
2478             }
2479         }
2480         GMX_LOG(mdlog.info).appendTextFormatted(
2481                 "PME domain decomposition: %d x %d x %d",
2482                 comm->npmenodes_x, comm->npmenodes_y, 1);
2483     }
2484     else
2485     {
2486         comm->npmedecompdim = 0;
2487         comm->npmenodes_x   = 0;
2488         comm->npmenodes_y   = 0;
2489     }
2490
2491     snew(comm->slb_frac, DIM);
2492     if (isDlbDisabled(comm))
2493     {
2494         comm->slb_frac[XX] = get_slb_frac(mdlog, "x", dd->nc[XX], options.cellSizeX);
2495         comm->slb_frac[YY] = get_slb_frac(mdlog, "y", dd->nc[YY], options.cellSizeY);
2496         comm->slb_frac[ZZ] = get_slb_frac(mdlog, "z", dd->nc[ZZ], options.cellSizeZ);
2497     }
2498
2499     /* Set the multi-body cut-off and cellsize limit for DLB */
2500     comm->cutoff_mbody   = systemInfo.minCutoffForMultiBody;
2501     comm->cellsize_limit = systemInfo.cellsizeLimit;
2502     if (systemInfo.haveInterDomainBondeds && systemInfo.increaseMultiBodyCutoff)
2503     {
2504         if (systemInfo.filterBondedCommunication || !isDlbDisabled(comm))
2505         {
2506             /* Set the bonded communication distance to halfway
2507              * the minimum and the maximum,
2508              * since the extra communication cost is nearly zero.
2509              */
2510             real acs           = average_cellsize_min(*ddbox, dd->nc);
2511             comm->cutoff_mbody = 0.5*(systemInfo.minCutoffForMultiBody + acs);
2512             if (!isDlbDisabled(comm))
2513             {
2514                 /* Check if this does not limit the scaling */
2515                 comm->cutoff_mbody = std::min(comm->cutoff_mbody,
2516                                               options.dlbScaling*acs);
2517             }
2518             if (!systemInfo.filterBondedCommunication)
2519             {
2520                 /* Without bBondComm do not go beyond the n.b. cut-off */
2521                 comm->cutoff_mbody = std::min(comm->cutoff_mbody, systemInfo.cutoff);
2522                 if (comm->cellsize_limit >= systemInfo.cutoff)
2523                 {
2524                     /* We don't loose a lot of efficieny
2525                      * when increasing it to the n.b. cut-off.
2526                      * It can even be slightly faster, because we need
2527                      * less checks for the communication setup.
2528                      */
2529                     comm->cutoff_mbody = systemInfo.cutoff;
2530                 }
2531             }
2532             /* Check if we did not end up below our original limit */
2533             comm->cutoff_mbody = std::max(comm->cutoff_mbody,
2534                                           systemInfo.minCutoffForMultiBody);
2535
2536             if (comm->cutoff_mbody > comm->cellsize_limit)
2537             {
2538                 comm->cellsize_limit = comm->cutoff_mbody;
2539             }
2540         }
2541         /* Without DLB and cutoff_mbody<cutoff, cutoff_mbody is dynamic */
2542     }
2543
2544     if (debug)
2545     {
2546         fprintf(debug, "Bonded atom communication beyond the cut-off: %s\n"
2547                 "cellsize limit %f\n",
2548                 gmx::boolToString(systemInfo.filterBondedCommunication),
2549                 comm->cellsize_limit);
2550     }
2551
2552     if (MASTER(cr))
2553     {
2554         check_dd_restrictions(dd, ir, mdlog);
2555     }
2556 }
2557
2558 static char *init_bLocalCG(const gmx_mtop_t *mtop)
2559 {
2560     int   ncg, cg;
2561     char *bLocalCG;
2562
2563     ncg = ncg_mtop(mtop);
2564     snew(bLocalCG, ncg);
2565     for (cg = 0; cg < ncg; cg++)
2566     {
2567         bLocalCG[cg] = FALSE;
2568     }
2569
2570     return bLocalCG;
2571 }
2572
2573 void dd_init_bondeds(FILE *fplog,
2574                      gmx_domdec_t *dd,
2575                      const gmx_mtop_t *mtop,
2576                      const gmx_vsite_t *vsite,
2577                      const t_inputrec *ir,
2578                      gmx_bool bBCheck, cginfo_mb_t *cginfo_mb)
2579 {
2580     gmx_domdec_comm_t *comm;
2581
2582     dd_make_reverse_top(fplog, dd, mtop, vsite, ir, bBCheck);
2583
2584     comm = dd->comm;
2585
2586     if (comm->systemInfo.filterBondedCommunication)
2587     {
2588         /* Communicate atoms beyond the cut-off for bonded interactions */
2589         comm = dd->comm;
2590
2591         comm->cglink = make_charge_group_links(mtop, dd, cginfo_mb);
2592
2593         comm->bLocalCG = init_bLocalCG(mtop);
2594     }
2595     else
2596     {
2597         /* Only communicate atoms based on cut-off */
2598         comm->cglink   = nullptr;
2599         comm->bLocalCG = nullptr;
2600     }
2601 }
2602
2603 static void writeSettings(gmx::TextWriter       *log,
2604                           gmx_domdec_t          *dd,
2605                           const gmx_mtop_t      *mtop,
2606                           const t_inputrec      *ir,
2607                           gmx_bool               bDynLoadBal,
2608                           real                   dlb_scale,
2609                           const gmx_ddbox_t     *ddbox)
2610 {
2611     gmx_domdec_comm_t *comm;
2612     int                d;
2613     ivec               np;
2614     real               limit, shrink;
2615
2616     comm = dd->comm;
2617
2618     if (bDynLoadBal)
2619     {
2620         log->writeString("The maximum number of communication pulses is:");
2621         for (d = 0; d < dd->ndim; d++)
2622         {
2623             log->writeStringFormatted(" %c %d", dim2char(dd->dim[d]), comm->cd[d].np_dlb);
2624         }
2625         log->ensureLineBreak();
2626         log->writeLineFormatted("The minimum size for domain decomposition cells is %.3f nm", comm->cellsize_limit);
2627         log->writeLineFormatted("The requested allowed shrink of DD cells (option -dds) is: %.2f", dlb_scale);
2628         log->writeString("The allowed shrink of domain decomposition cells is:");
2629         for (d = 0; d < DIM; d++)
2630         {
2631             if (dd->nc[d] > 1)
2632             {
2633                 if (d >= ddbox->npbcdim && dd->nc[d] == 2)
2634                 {
2635                     shrink = 0;
2636                 }
2637                 else
2638                 {
2639                     shrink =
2640                         comm->cellsize_min_dlb[d]/
2641                         (ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
2642                 }
2643                 log->writeStringFormatted(" %c %.2f", dim2char(d), shrink);
2644             }
2645         }
2646         log->ensureLineBreak();
2647     }
2648     else
2649     {
2650         set_dd_cell_sizes_slb(dd, ddbox, setcellsizeslbPULSE_ONLY, np);
2651         log->writeString("The initial number of communication pulses is:");
2652         for (d = 0; d < dd->ndim; d++)
2653         {
2654             log->writeStringFormatted(" %c %d", dim2char(dd->dim[d]), np[dd->dim[d]]);
2655         }
2656         log->ensureLineBreak();
2657         log->writeString("The initial domain decomposition cell size is:");
2658         for (d = 0; d < DIM; d++)
2659         {
2660             if (dd->nc[d] > 1)
2661             {
2662                 log->writeStringFormatted(" %c %.2f nm",
2663                                           dim2char(d), dd->comm->cellsize_min[d]);
2664             }
2665         }
2666         log->ensureLineBreak();
2667         log->writeLine();
2668     }
2669
2670     const bool haveInterDomainVsites =
2671         (countInterUpdategroupVsites(*mtop, comm->systemInfo.updateGroupingPerMoleculetype) != 0);
2672
2673     if (comm->systemInfo.haveInterDomainBondeds ||
2674         haveInterDomainVsites ||
2675         comm->systemInfo.haveSplitConstraints ||
2676         comm->systemInfo.haveSplitSettles)
2677     {
2678         std::string decompUnits;
2679         if (comm->systemInfo.useUpdateGroups)
2680         {
2681             decompUnits = "atom groups";
2682         }
2683         else
2684         {
2685             decompUnits = "atoms";
2686         }
2687
2688         log->writeLineFormatted("The maximum allowed distance for %s involved in interactions is:", decompUnits.c_str());
2689         log->writeLineFormatted("%40s  %-7s %6.3f nm", "non-bonded interactions", "", comm->systemInfo.cutoff);
2690
2691         if (bDynLoadBal)
2692         {
2693             limit = dd->comm->cellsize_limit;
2694         }
2695         else
2696         {
2697             if (dd->unitCellInfo.ddBoxIsDynamic)
2698             {
2699                 log->writeLine("(the following are initial values, they could change due to box deformation)");
2700             }
2701             limit = dd->comm->cellsize_min[XX];
2702             for (d = 1; d < DIM; d++)
2703             {
2704                 limit = std::min(limit, dd->comm->cellsize_min[d]);
2705             }
2706         }
2707
2708         if (comm->systemInfo.haveInterDomainBondeds)
2709         {
2710             log->writeLineFormatted("%40s  %-7s %6.3f nm",
2711                                     "two-body bonded interactions", "(-rdd)",
2712                                     std::max(comm->systemInfo.cutoff, comm->cutoff_mbody));
2713             log->writeLineFormatted("%40s  %-7s %6.3f nm",
2714                                     "multi-body bonded interactions", "(-rdd)",
2715                                     (comm->systemInfo.filterBondedCommunication || isDlbOn(dd->comm)) ? comm->cutoff_mbody : std::min(comm->systemInfo.cutoff, limit));
2716         }
2717         if (haveInterDomainVsites)
2718         {
2719             log->writeLineFormatted("%40s  %-7s %6.3f nm",
2720                                     "virtual site constructions", "(-rcon)", limit);
2721         }
2722         if (comm->systemInfo.haveSplitConstraints || comm->systemInfo.haveSplitSettles)
2723         {
2724             std::string separation = gmx::formatString("atoms separated by up to %d constraints",
2725                                                        1+ir->nProjOrder);
2726             log->writeLineFormatted("%40s  %-7s %6.3f nm\n",
2727                                     separation.c_str(), "(-rcon)", limit);
2728         }
2729         log->ensureLineBreak();
2730     }
2731 }
2732
2733 static void logSettings(const gmx::MDLogger &mdlog,
2734                         gmx_domdec_t        *dd,
2735                         const gmx_mtop_t    *mtop,
2736                         const t_inputrec    *ir,
2737                         real                 dlb_scale,
2738                         const gmx_ddbox_t   *ddbox)
2739 {
2740     gmx::StringOutputStream stream;
2741     gmx::TextWriter         log(&stream);
2742     writeSettings(&log, dd, mtop, ir, isDlbOn(dd->comm), dlb_scale, ddbox);
2743     if (dd->comm->dlbState == DlbState::offCanTurnOn)
2744     {
2745         {
2746             log.ensureEmptyLine();
2747             log.writeLine("When dynamic load balancing gets turned on, these settings will change to:");
2748         }
2749         writeSettings(&log, dd, mtop, ir, true, dlb_scale, ddbox);
2750     }
2751     GMX_LOG(mdlog.info).asParagraph().appendText(stream.toString());
2752 }
2753
2754 static void set_cell_limits_dlb(const gmx::MDLogger &mdlog,
2755                                 gmx_domdec_t        *dd,
2756                                 real                 dlb_scale,
2757                                 const t_inputrec    *ir,
2758                                 const gmx_ddbox_t   *ddbox)
2759 {
2760     gmx_domdec_comm_t *comm;
2761     int                d, dim, npulse, npulse_d_max, npulse_d;
2762     gmx_bool           bNoCutOff;
2763
2764     comm = dd->comm;
2765
2766     bNoCutOff = (ir->rvdw == 0 || ir->rcoulomb == 0);
2767
2768     /* Determine the maximum number of comm. pulses in one dimension */
2769
2770     comm->cellsize_limit = std::max(comm->cellsize_limit, comm->cutoff_mbody);
2771
2772     /* Determine the maximum required number of grid pulses */
2773     if (comm->cellsize_limit >= comm->systemInfo.cutoff)
2774     {
2775         /* Only a single pulse is required */
2776         npulse = 1;
2777     }
2778     else if (!bNoCutOff && comm->cellsize_limit > 0)
2779     {
2780         /* We round down slightly here to avoid overhead due to the latency
2781          * of extra communication calls when the cut-off
2782          * would be only slightly longer than the cell size.
2783          * Later cellsize_limit is redetermined,
2784          * so we can not miss interactions due to this rounding.
2785          */
2786         npulse = static_cast<int>(0.96 + comm->systemInfo.cutoff/comm->cellsize_limit);
2787     }
2788     else
2789     {
2790         /* There is no cell size limit */
2791         npulse = std::max(dd->nc[XX]-1, std::max(dd->nc[YY]-1, dd->nc[ZZ]-1));
2792     }
2793
2794     if (!bNoCutOff && npulse > 1)
2795     {
2796         /* See if we can do with less pulses, based on dlb_scale */
2797         npulse_d_max = 0;
2798         for (d = 0; d < dd->ndim; d++)
2799         {
2800             dim      = dd->dim[d];
2801             npulse_d = static_cast<int>(1 + dd->nc[dim]*comm->systemInfo.cutoff
2802                                         /(ddbox->box_size[dim]*ddbox->skew_fac[dim]*dlb_scale));
2803             npulse_d_max = std::max(npulse_d_max, npulse_d);
2804         }
2805         npulse = std::min(npulse, npulse_d_max);
2806     }
2807
2808     /* This env var can override npulse */
2809     d = dd_getenv(mdlog, "GMX_DD_NPULSE", 0);
2810     if (d > 0)
2811     {
2812         npulse = d;
2813     }
2814
2815     comm->maxpulse       = 1;
2816     comm->bVacDLBNoLimit = (ir->ePBC == epbcNONE);
2817     for (d = 0; d < dd->ndim; d++)
2818     {
2819         comm->cd[d].np_dlb    = std::min(npulse, dd->nc[dd->dim[d]]-1);
2820         comm->maxpulse        = std::max(comm->maxpulse, comm->cd[d].np_dlb);
2821         if (comm->cd[d].np_dlb < dd->nc[dd->dim[d]]-1)
2822         {
2823             comm->bVacDLBNoLimit = FALSE;
2824         }
2825     }
2826
2827     /* cellsize_limit is set for LINCS in init_domain_decomposition */
2828     if (!comm->bVacDLBNoLimit)
2829     {
2830         comm->cellsize_limit = std::max(comm->cellsize_limit,
2831                                         comm->systemInfo.cutoff/comm->maxpulse);
2832     }
2833     comm->cellsize_limit = std::max(comm->cellsize_limit, comm->cutoff_mbody);
2834     /* Set the minimum cell size for each DD dimension */
2835     for (d = 0; d < dd->ndim; d++)
2836     {
2837         if (comm->bVacDLBNoLimit ||
2838             comm->cd[d].np_dlb*comm->cellsize_limit >= comm->systemInfo.cutoff)
2839         {
2840             comm->cellsize_min_dlb[dd->dim[d]] = comm->cellsize_limit;
2841         }
2842         else
2843         {
2844             comm->cellsize_min_dlb[dd->dim[d]] =
2845                 comm->systemInfo.cutoff/comm->cd[d].np_dlb;
2846         }
2847     }
2848     if (comm->cutoff_mbody <= 0)
2849     {
2850         comm->cutoff_mbody = std::min(comm->systemInfo.cutoff, comm->cellsize_limit);
2851     }
2852     if (isDlbOn(comm))
2853     {
2854         set_dlb_limits(dd);
2855     }
2856 }
2857
2858 gmx_bool dd_bonded_molpbc(const gmx_domdec_t *dd, int ePBC)
2859 {
2860     /* If each molecule is a single charge group
2861      * or we use domain decomposition for each periodic dimension,
2862      * we do not need to take pbc into account for the bonded interactions.
2863      */
2864     return (ePBC != epbcNONE && dd->comm->systemInfo.haveInterDomainBondeds &&
2865             !(dd->nc[XX] > 1 &&
2866               dd->nc[YY] > 1 &&
2867               (dd->nc[ZZ] > 1 || ePBC == epbcXY)));
2868 }
2869
2870 /*! \brief Sets grid size limits and PP-PME setup, prints settings to log */
2871 static void set_ddgrid_parameters(const gmx::MDLogger &mdlog,
2872                                   gmx_domdec_t *dd, real dlb_scale,
2873                                   const gmx_mtop_t *mtop, const t_inputrec *ir,
2874                                   const gmx_ddbox_t *ddbox)
2875 {
2876     gmx_domdec_comm_t *comm;
2877     int                natoms_tot;
2878     real               vol_frac;
2879
2880     comm = dd->comm;
2881
2882     if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
2883     {
2884         init_ddpme(dd, &comm->ddpme[0], 0);
2885         if (comm->npmedecompdim >= 2)
2886         {
2887             init_ddpme(dd, &comm->ddpme[1], 1);
2888         }
2889     }
2890     else
2891     {
2892         comm->npmenodes = 0;
2893         if (dd->pme_nodeid >= 0)
2894         {
2895             gmx_fatal_collective(FARGS, dd->mpi_comm_all, DDMASTER(dd),
2896                                  "Can not have separate PME ranks without PME electrostatics");
2897         }
2898     }
2899
2900     if (debug)
2901     {
2902         fprintf(debug, "The DD cut-off is %f\n", comm->systemInfo.cutoff);
2903     }
2904     if (!isDlbDisabled(comm))
2905     {
2906         set_cell_limits_dlb(mdlog, dd, dlb_scale, ir, ddbox);
2907     }
2908
2909     logSettings(mdlog, dd, mtop, ir, dlb_scale, ddbox);
2910
2911     if (ir->ePBC == epbcNONE)
2912     {
2913         vol_frac = 1 - 1/static_cast<double>(dd->nnodes);
2914     }
2915     else
2916     {
2917         vol_frac =
2918             (1 + comm_box_frac(dd->nc, comm->systemInfo.cutoff, ddbox))/static_cast<double>(dd->nnodes);
2919     }
2920     if (debug)
2921     {
2922         fprintf(debug, "Volume fraction for all DD zones: %f\n", vol_frac);
2923     }
2924     natoms_tot = comm->cgs_gl.index[comm->cgs_gl.nr];
2925
2926     dd->ga2la  = new gmx_ga2la_t(natoms_tot,
2927                                  static_cast<int>(vol_frac*natoms_tot));
2928 }
2929
2930 /*! \brief Get some important DD parameters which can be modified by env.vars */
2931 static DDSettings
2932 getDDSettings(const gmx::MDLogger     &mdlog,
2933               const DomdecOptions     &options,
2934               const gmx::MdrunOptions &mdrunOptions,
2935               const t_inputrec        &ir)
2936 {
2937     DDSettings ddSettings;
2938
2939     ddSettings.useSendRecv2  = (dd_getenv(mdlog, "GMX_DD_USE_SENDRECV2", 0) != 0);
2940     ddSettings.dlb_scale_lim = dd_getenv(mdlog, "GMX_DLB_MAX_BOX_SCALING", 10);
2941     ddSettings.eFlop         = dd_getenv(mdlog, "GMX_DLB_BASED_ON_FLOPS", 0);
2942     const int recload        = dd_getenv(mdlog, "GMX_DD_RECORD_LOAD", 1);
2943     ddSettings.nstDDDump     = dd_getenv(mdlog, "GMX_DD_NST_DUMP", 0);
2944     ddSettings.nstDDDumpGrid = dd_getenv(mdlog, "GMX_DD_NST_DUMP_GRID", 0);
2945     ddSettings.DD_debug      = dd_getenv(mdlog, "GMX_DD_DEBUG", 0);
2946
2947     if (ddSettings.useSendRecv2)
2948     {
2949         GMX_LOG(mdlog.info).appendText("Will use two sequential MPI_Sendrecv calls instead of two simultaneous non-blocking MPI_Irecv and MPI_Isend pairs for constraint and vsite communication");
2950     }
2951
2952     if (ddSettings.eFlop)
2953     {
2954         GMX_LOG(mdlog.info).appendText("Will load balance based on FLOP count");
2955         ddSettings.recordLoad = true;
2956     }
2957     else
2958     {
2959         ddSettings.recordLoad = (wallcycle_have_counter() && recload > 0);
2960     }
2961
2962     ddSettings.initialDlbState =
2963         determineInitialDlbState(mdlog, options.dlbOption, ddSettings.recordLoad, mdrunOptions, &ir);
2964
2965     return ddSettings;
2966 }
2967
2968 gmx_domdec_t::gmx_domdec_t(const t_inputrec &ir) :
2969     unitCellInfo(ir)
2970 {
2971 }
2972
2973 gmx_domdec_t *init_domain_decomposition(const gmx::MDLogger           &mdlog,
2974                                         t_commrec                     *cr,
2975                                         const DomdecOptions           &options,
2976                                         const gmx::MdrunOptions       &mdrunOptions,
2977                                         const gmx_mtop_t              *mtop,
2978                                         const t_inputrec              *ir,
2979                                         const matrix                   box,
2980                                         gmx::ArrayRef<const gmx::RVec> xGlobal,
2981                                         gmx::LocalAtomSetManager      *atomSets)
2982 {
2983     gmx_domdec_t      *dd;
2984
2985     GMX_LOG(mdlog.info).appendTextFormatted(
2986             "\nInitializing Domain Decomposition on %d ranks", cr->nnodes);
2987
2988     dd = new gmx_domdec_t(*ir);
2989
2990     dd->comm = init_dd_comm();
2991
2992     DDSettings  ddSettings = getDDSettings(mdlog, options, mdrunOptions, *ir);
2993     if (ddSettings.eFlop > 1)
2994     {
2995         /* Ensure that we have different random flop counts on different ranks */
2996         srand(1 + cr->nodeid);
2997     }
2998
2999     gmx_ddbox_t ddbox = {0};
3000     set_dd_limits_and_grid(mdlog, cr, dd, options,
3001                            ddSettings,
3002                            mtop, ir,
3003                            box, xGlobal,
3004                            &ddbox);
3005
3006     make_dd_communicators(mdlog, cr, dd, options.rankOrder);
3007
3008     if (thisRankHasDuty(cr, DUTY_PP))
3009     {
3010         set_ddgrid_parameters(mdlog, dd, options.dlbScaling, mtop, ir, &ddbox);
3011
3012         setup_neighbor_relations(dd);
3013     }
3014
3015     /* Set overallocation to avoid frequent reallocation of arrays */
3016     set_over_alloc_dd(TRUE);
3017
3018     dd->atomSets = atomSets;
3019
3020     return dd;
3021 }
3022
3023 static gmx_bool test_dd_cutoff(t_commrec     *cr,
3024                                const t_state &state,
3025                                real           cutoffRequested)
3026 {
3027     gmx_domdec_t *dd;
3028     gmx_ddbox_t   ddbox;
3029     int           d, dim, np;
3030     real          inv_cell_size;
3031     int           LocallyLimited;
3032
3033     dd = cr->dd;
3034
3035     set_ddbox(*dd, false, state.box, true, state.x, &ddbox);
3036
3037     LocallyLimited = 0;
3038
3039     for (d = 0; d < dd->ndim; d++)
3040     {
3041         dim = dd->dim[d];
3042
3043         inv_cell_size = DD_CELL_MARGIN*dd->nc[dim]/ddbox.box_size[dim];
3044         if (dd->unitCellInfo.ddBoxIsDynamic)
3045         {
3046             inv_cell_size *= DD_PRES_SCALE_MARGIN;
3047         }
3048
3049         np = 1 + static_cast<int>(cutoffRequested*inv_cell_size*ddbox.skew_fac[dim]);
3050
3051         if (!isDlbDisabled(dd->comm) && (dim < ddbox.npbcdim) && (dd->comm->cd[d].np_dlb > 0))
3052         {
3053             if (np > dd->comm->cd[d].np_dlb)
3054             {
3055                 return FALSE;
3056             }
3057
3058             /* If a current local cell size is smaller than the requested
3059              * cut-off, we could still fix it, but this gets very complicated.
3060              * Without fixing here, we might actually need more checks.
3061              */
3062             real cellSizeAlongDim = (dd->comm->cell_x1[dim] - dd->comm->cell_x0[dim])*ddbox.skew_fac[dim];
3063             if (cellSizeAlongDim*dd->comm->cd[d].np_dlb < cutoffRequested)
3064             {
3065                 LocallyLimited = 1;
3066             }
3067         }
3068     }
3069
3070     if (!isDlbDisabled(dd->comm))
3071     {
3072         /* If DLB is not active yet, we don't need to check the grid jumps.
3073          * Actually we shouldn't, because then the grid jump data is not set.
3074          */
3075         if (isDlbOn(dd->comm) &&
3076             check_grid_jump(0, dd, cutoffRequested, &ddbox, FALSE))
3077         {
3078             LocallyLimited = 1;
3079         }
3080
3081         gmx_sumi(1, &LocallyLimited, cr);
3082
3083         if (LocallyLimited > 0)
3084         {
3085             return FALSE;
3086         }
3087     }
3088
3089     return TRUE;
3090 }
3091
3092 gmx_bool change_dd_cutoff(t_commrec     *cr,
3093                           const t_state &state,
3094                           real           cutoffRequested)
3095 {
3096     gmx_bool bCutoffAllowed;
3097
3098     bCutoffAllowed = test_dd_cutoff(cr, state, cutoffRequested);
3099
3100     if (bCutoffAllowed)
3101     {
3102         cr->dd->comm->systemInfo.cutoff = cutoffRequested;
3103     }
3104
3105     return bCutoffAllowed;
3106 }