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                ddnodeid = -1;
 165
 166     const DDRankSetup &ddRankSetup = dd->comm->ddRankSetup;
 167     const int          ddindex     = dd_index(dd->nc, c);
 168     if (ddRankSetup.bCartesianPP_PME)
 169     {
 170         ddnodeid = ddRankSetup.ddindex2ddnodeid[ddindex];
 171     }
 172     else if (ddRankSetup.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 DDRankSetup &ddRankSetup,
 695                                    const ivec         coord,
 696                                    ivec               coord_pme)
 697 {
 698     const int nc   = ddRankSetup.numPPCells[ddRankSetup.cartpmedim];
 699     const int ntot = ddRankSetup.ntot[ddRankSetup.cartpmedim];
 700     copy_ivec(coord, coord_pme);
 701     coord_pme[ddRankSetup.cartpmedim] =
 702         nc + (coord[ddRankSetup.cartpmedim]*(ntot - nc) + (ntot - nc)/2)/nc;
 703 }
 704
 705 /* Returns the PME rank index in 0...npmenodes-1 for the PP cell with index ddCellIndex */
 706 static int ddindex2pmeindex(const DDRankSetup &ddRankSetup,
 707                             const int          ddCellIndex)
 708 {
 709     const int npp  = ddRankSetup.numPPRanks;
 710     const int npme = ddRankSetup.npmenodes;
 711
 712     /* Here we assign a PME node to communicate with this DD node
 713      * by assuming that the major index of both is x.
 714      * We add npme/2 to obtain an even distribution.
 715      */
 716     return (ddCellIndex*npme + npme/2)/npp;
 717 }
 718
 719 static int *dd_interleaved_pme_ranks(const DDRankSetup &ddRankSetup)
 720 {
 721     int *pme_rank;
 722     int  n, i, p0, p1;
 723
 724     snew(pme_rank, ddRankSetup.npmenodes);
 725     n = 0;
 726     for (i = 0; i < ddRankSetup.numPPRanks; i++)
 727     {
 728         p0 = ddindex2pmeindex(ddRankSetup, i);
 729         p1 = ddindex2pmeindex(ddRankSetup, i + 1);
 730         if (i + 1 == ddRankSetup.numPPRanks || p1 > p0)
 731         {
 732             if (debug)
 733             {
 734                 fprintf(debug, "pme_rank[%d] = %d\n", n, i+1+n);
 735             }
 736             pme_rank[n] = i + 1 + n;
 737             n++;
 738         }
 739     }
 740
 741     return pme_rank;
 742 }
 743
 744 static int gmx_ddcoord2pmeindex(const t_commrec *cr, int x, int y, int z)
 745 {
 746     gmx_domdec_t *dd;
 747     ivec          coords;
 748     int           slab;
 749
 750     dd = cr->dd;
 751     /*
 752        if (dd->comm->bCartesian) {
 753        gmx_ddindex2xyz(dd->nc,ddindex,coords);
 754        dd_coords2pmecoords(dd,coords,coords_pme);
 755        copy_ivec(dd->ntot,nc);
 756        nc[dd->cartpmedim]         -= dd->nc[dd->cartpmedim];
 757        coords_pme[dd->cartpmedim] -= dd->nc[dd->cartpmedim];
 758
 759        slab = (coords_pme[XX]*nc[YY] + coords_pme[YY])*nc[ZZ] + coords_pme[ZZ];
 760        } else {
 761        slab = (ddindex*cr->npmenodes + cr->npmenodes/2)/dd->nnodes;
 762        }
 763      */
 764     coords[XX] = x;
 765     coords[YY] = y;
 766     coords[ZZ] = z;
 767     slab       = ddindex2pmeindex(dd->comm->ddRankSetup, dd_index(dd->nc, coords));
 768
 769     return slab;
 770 }
 771
 772 static int ddcoord2simnodeid(const t_commrec *cr, int x, int y, int z)
 773 {
 774     const DDRankSetup &ddRankSetup = cr->dd->comm->ddRankSetup;
 775     ivec               coords;
 776     int                ddindex, nodeid = -1;
 777
 778     coords[XX] = x;
 779     coords[YY] = y;
 780     coords[ZZ] = z;
 781     if (ddRankSetup.bCartesianPP_PME)
 782     {
 783 #if GMX_MPI
 784         MPI_Cart_rank(cr->mpi_comm_mysim, coords, &nodeid);
 785 #endif
 786     }
 787     else
 788     {
 789         ddindex = dd_index(cr->dd->nc, coords);
 790         if (ddRankSetup.bCartesianPP)
 791         {
 792             nodeid = ddRankSetup.ddindex2simnodeid[ddindex];
 793         }
 794         else
 795         {
 796             if (ddRankSetup.pmenodes)
 797             {
 798                 nodeid = ddindex + gmx_ddcoord2pmeindex(cr, x, y, z);
 799             }
 800             else
 801             {
 802                 nodeid = ddindex;
 803             }
 804         }
 805     }
 806
 807     return nodeid;
 808 }
 809
 810 static int dd_simnode2pmenode(const DDRankSetup          &ddRankSetup,
 811                               const t_commrec gmx_unused *cr,
 812                               const int                   sim_nodeid)
 813 {
 814     int                pmenode = -1;
 815
 816     /* This assumes a uniform x domain decomposition grid cell size */
 817     if (ddRankSetup.bCartesianPP_PME)
 818     {
 819 #if GMX_MPI
 820         ivec coord, coord_pme;
 821         MPI_Cart_coords(cr->mpi_comm_mysim, sim_nodeid, DIM, coord);
 822         if (coord[ddRankSetup.cartpmedim] < ddRankSetup.numPPCells[ddRankSetup.cartpmedim])
 823         {
 824             /* This is a PP rank */
 825             dd_cart_coord2pmecoord(ddRankSetup, coord, coord_pme);
 826             MPI_Cart_rank(cr->mpi_comm_mysim, coord_pme, &pmenode);
 827         }
 828 #endif
 829     }
 830     else if (ddRankSetup.bCartesianPP)
 831     {
 832         if (sim_nodeid < ddRankSetup.numPPRanks)
 833         {
 834             pmenode = ddRankSetup.numPPRanks + ddindex2pmeindex(ddRankSetup, sim_nodeid);
 835         }
 836     }
 837     else
 838     {
 839         /* This assumes DD cells with identical x coordinates
 840          * are numbered sequentially.
 841          */
 842         if (ddRankSetup.pmenodes == nullptr)
 843         {
 844             if (sim_nodeid < ddRankSetup.numPPRanks)
 845             {
 846                 /* The DD index equals the nodeid */
 847                 pmenode = ddRankSetup.numPPRanks + ddindex2pmeindex(ddRankSetup, sim_nodeid);
 848             }
 849         }
 850         else
 851         {
 852             int i = 0;
 853             while (sim_nodeid > ddRankSetup.pmenodes[i])
 854             {
 855                 i++;
 856             }
 857             if (sim_nodeid < ddRankSetup.pmenodes[i])
 858             {
 859                 pmenode = ddRankSetup.pmenodes[i];
 860             }
 861         }
 862     }
 863
 864     return pmenode;
 865 }
 866
 867 NumPmeDomains getNumPmeDomains(const gmx_domdec_t *dd)
 868 {
 869     if (dd != nullptr)
 870     {
 871         return {
 872                    dd->comm->ddRankSetup.npmenodes_x,
 873                    dd->comm->ddRankSetup.npmenodes_y
 874         };
 875     }
 876     else
 877     {
 878         return {
 879                    1, 1
 880         };
 881     }
 882 }
 883
 884 std::vector<int> get_pme_ddranks(const t_commrec *cr,
 885                                  const int        pmenodeid)
 886 {
 887     const DDRankSetup &ddRankSetup = cr->dd->comm->ddRankSetup;
 888     GMX_RELEASE_ASSERT(ddRankSetup.npmenodes > 0, "This function should only be called when PME-only ranks are in use");
 889     std::vector<int>   ddranks;
 890     ddranks.reserve((ddRankSetup.numPPRanks + ddRankSetup.npmenodes - 1)/ddRankSetup.npmenodes);
 891
 892     for (int x = 0; x < ddRankSetup.numPPCells[XX]; x++)
 893     {
 894         for (int y = 0; y < ddRankSetup.numPPCells[YY]; y++)
 895         {
 896             for (int z = 0; z < ddRankSetup.numPPCells[ZZ]; z++)
 897             {
 898                 if (ddRankSetup.bCartesianPP_PME)
 899                 {
 900                     ivec coord = { x, y, z };
 901                     ivec coord_pme;
 902                     dd_cart_coord2pmecoord(ddRankSetup, coord, coord_pme);
 903                     if (cr->dd->ci[XX] == coord_pme[XX] &&
 904                         cr->dd->ci[YY] == coord_pme[YY] &&
 905                         cr->dd->ci[ZZ] == coord_pme[ZZ])
 906                     {
 907                         ddranks.push_back(ddcoord2simnodeid(cr, x, y, z));
 908                     }
 909                 }
 910                 else
 911                 {
 912                     /* The slab corresponds to the nodeid in the PME group */
 913                     if (gmx_ddcoord2pmeindex(cr, x, y, z) == pmenodeid)
 914                     {
 915                         ddranks.push_back(ddcoord2simnodeid(cr, x, y, z));
 916                     }
 917                 }
 918             }
 919         }
 920     }
 921     return ddranks;
 922 }
 923
 924 static gmx_bool receive_vir_ener(const gmx_domdec_t *dd, const t_commrec *cr)
 925 {
 926     gmx_bool bReceive = TRUE;
 927
 928     if (cr->npmenodes < dd->nnodes)
 929     {
 930         const DDRankSetup &ddRankSetup = dd->comm->ddRankSetup;
 931         if (ddRankSetup.bCartesianPP_PME)
 932         {
 933 #if GMX_MPI
 934             int  pmenode = dd_simnode2pmenode(ddRankSetup, cr, cr->sim_nodeid);
 935             ivec coords;
 936             MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, coords);
 937             coords[ddRankSetup.cartpmedim]++;
 938             if (coords[ddRankSetup.cartpmedim] < dd->nc[ddRankSetup.cartpmedim])
 939             {
 940                 int rank;
 941                 MPI_Cart_rank(cr->mpi_comm_mysim, coords, &rank);
 942                 if (dd_simnode2pmenode(ddRankSetup, cr, rank) == pmenode)
 943                 {
 944                     /* This is not the last PP node for pmenode */
 945                     bReceive = FALSE;
 946                 }
 947             }
 948 #else
 949             GMX_RELEASE_ASSERT(false, "Without MPI we should not have Cartesian PP-PME with #PMEnodes < #DDnodes");
 950 #endif
 951         }
 952         else
 953         {
 954             int pmenode = dd_simnode2pmenode(ddRankSetup, cr, cr->sim_nodeid);
 955             if (cr->sim_nodeid+1 < cr->nnodes &&
 956                 dd_simnode2pmenode(ddRankSetup, cr, cr->sim_nodeid + 1) == pmenode)
 957             {
 958                 /* This is not the last PP node for pmenode */
 959                 bReceive = FALSE;
 960             }
 961         }
 962     }
 963
 964     return bReceive;
 965 }
 966
 967 static void set_slb_pme_dim_f(gmx_domdec_t *dd, int dim, real **dim_f)
 968 {
 969     gmx_domdec_comm_t *comm;
 970     int                i;
 971
 972     comm = dd->comm;
 973
 974     snew(*dim_f, dd->nc[dim]+1);
 975     (*dim_f)[0] = 0;
 976     for (i = 1; i < dd->nc[dim]; i++)
 977     {
 978         if (comm->slb_frac[dim])
 979         {
 980             (*dim_f)[i] = (*dim_f)[i-1] + comm->slb_frac[dim][i-1];
 981         }
 982         else
 983         {
 984             (*dim_f)[i] = static_cast<real>(i)/static_cast<real>(dd->nc[dim]);
 985         }
 986     }
 987     (*dim_f)[dd->nc[dim]] = 1;
 988 }
 989
 990 static void init_ddpme(gmx_domdec_t *dd, gmx_ddpme_t *ddpme, int dimind)
 991 {
 992     const DDRankSetup &ddRankSetup = dd->comm->ddRankSetup;
 993
 994     if (dimind == 0 && dd->dim[0] == YY && ddRankSetup.npmenodes_x == 1)
 995     {
 996         ddpme->dim = YY;
 997     }
 998     else
 999     {
1000         ddpme->dim = dimind;
1001     }
1002     ddpme->dim_match = (ddpme->dim == dd->dim[dimind]);
1003
1004     ddpme->nslab = (ddpme->dim == 0 ?
1005                     ddRankSetup.npmenodes_x :
1006                     ddRankSetup.npmenodes_y);
1007
1008     if (ddpme->nslab <= 1)
1009     {
1010         return;
1011     }
1012
1013     const int nso = ddRankSetup.npmenodes/ddpme->nslab;
1014     /* Determine for each PME slab the PP location range for dimension dim */
1015     snew(ddpme->pp_min, ddpme->nslab);
1016     snew(ddpme->pp_max, ddpme->nslab);
1017     for (int slab = 0; slab < ddpme->nslab; slab++)
1018     {
1019         ddpme->pp_min[slab] = dd->nc[dd->dim[dimind]] - 1;
1020         ddpme->pp_max[slab] = 0;
1021     }
1022     for (int i = 0; i < dd->nnodes; i++)
1023     {
1024         ivec xyz;
1025         ddindex2xyz(dd->nc, i, xyz);
1026         /* For y only use our y/z slab.
1027          * This assumes that the PME x grid size matches the DD grid size.
1028          */
1029         if (dimind == 0 || xyz[XX] == dd->ci[XX])
1030         {
1031             const int pmeindex = ddindex2pmeindex(ddRankSetup, i);
1032             int       slab;
1033             if (dimind == 0)
1034             {
1035                 slab = pmeindex/nso;
1036             }
1037             else
1038             {
1039                 slab = pmeindex % ddpme->nslab;
1040             }
1041             ddpme->pp_min[slab] = std::min(ddpme->pp_min[slab], xyz[dimind]);
1042             ddpme->pp_max[slab] = std::max(ddpme->pp_max[slab], xyz[dimind]);
1043         }
1044     }
1045
1046     set_slb_pme_dim_f(dd, ddpme->dim, &ddpme->slb_dim_f);
1047 }
1048
1049 int dd_pme_maxshift_x(const gmx_domdec_t *dd)
1050 {
1051     const DDRankSetup &ddRankSetup = dd->comm->ddRankSetup;
1052
1053     if (ddRankSetup.ddpme[0].dim == XX)
1054     {
1055         return ddRankSetup.ddpme[0].maxshift;
1056     }
1057     else
1058     {
1059         return 0;
1060     }
1061 }
1062
1063 int dd_pme_maxshift_y(const gmx_domdec_t *dd)
1064 {
1065     const DDRankSetup &ddRankSetup = dd->comm->ddRankSetup;
1066
1067     if (ddRankSetup.ddpme[0].dim == YY)
1068     {
1069         return ddRankSetup.ddpme[0].maxshift;
1070     }
1071     else if (ddRankSetup.npmedecompdim >= 2 && ddRankSetup.ddpme[1].dim == YY)
1072     {
1073         return ddRankSetup.ddpme[1].maxshift;
1074     }
1075     else
1076     {
1077         return 0;
1078     }
1079 }
1080
1081 bool ddHaveSplitConstraints(const gmx_domdec_t &dd)
1082 {
1083     return dd.comm->systemInfo.haveSplitConstraints;
1084 }
1085
1086 void dd_cycles_add(const gmx_domdec_t *dd, float cycles, int ddCycl)
1087 {
1088     /* Note that the cycles value can be incorrect, either 0 or some
1089      * extremely large value, when our thread migrated to another core
1090      * with an unsynchronized cycle counter. If this happens less often
1091      * that once per nstlist steps, this will not cause issues, since
1092      * we later subtract the maximum value from the sum over nstlist steps.
1093      * A zero count will slightly lower the total, but that's a small effect.
1094      * Note that the main purpose of the subtraction of the maximum value
1095      * is to avoid throwing off the load balancing when stalls occur due
1096      * e.g. system activity or network congestion.
1097      */
1098     dd->comm->cycl[ddCycl] += cycles;
1099     dd->comm->cycl_n[ddCycl]++;
1100     if (cycles > dd->comm->cycl_max[ddCycl])
1101     {
1102         dd->comm->cycl_max[ddCycl] = cycles;
1103     }
1104 }
1105
1106 #if GMX_MPI
1107 static void make_load_communicator(gmx_domdec_t *dd, int dim_ind, ivec loc)
1108 {
1109     MPI_Comm           c_row;
1110     int                dim, i, rank;
1111     ivec               loc_c;
1112     gmx_bool           bPartOfGroup = FALSE;
1113
1114     dim = dd->dim[dim_ind];
1115     copy_ivec(loc, loc_c);
1116     for (i = 0; i < dd->nc[dim]; i++)
1117     {
1118         loc_c[dim] = i;
1119         rank       = dd_index(dd->nc, loc_c);
1120         if (rank == dd->rank)
1121         {
1122             /* This process is part of the group */
1123             bPartOfGroup = TRUE;
1124         }
1125     }
1126     MPI_Comm_split(dd->mpi_comm_all, bPartOfGroup ? 0 : MPI_UNDEFINED, dd->rank,
1127                    &c_row);
1128     if (bPartOfGroup)
1129     {
1130         dd->comm->mpi_comm_load[dim_ind] = c_row;
1131         if (!isDlbDisabled(dd->comm))
1132         {
1133             DDCellsizesWithDlb &cellsizes = dd->comm->cellsizesWithDlb[dim_ind];
1134
1135             if (dd->ci[dim] == dd->master_ci[dim])
1136             {
1137                 /* This is the root process of this row */
1138                 cellsizes.rowMaster  = std::make_unique<RowMaster>();
1139
1140                 RowMaster &rowMaster = *cellsizes.rowMaster;
1141                 rowMaster.cellFrac.resize(ddCellFractionBufferSize(dd, dim_ind));
1142                 rowMaster.oldCellFrac.resize(dd->nc[dim] + 1);
1143                 rowMaster.isCellMin.resize(dd->nc[dim]);
1144                 if (dim_ind > 0)
1145                 {
1146                     rowMaster.bounds.resize(dd->nc[dim]);
1147                 }
1148                 rowMaster.buf_ncd.resize(dd->nc[dim]);
1149             }
1150             else
1151             {
1152                 /* This is not a root process, we only need to receive cell_f */
1153                 cellsizes.fracRow.resize(ddCellFractionBufferSize(dd, dim_ind));
1154             }
1155         }
1156         if (dd->ci[dim] == dd->master_ci[dim])
1157         {
1158             snew(dd->comm->load[dim_ind].load, dd->nc[dim]*DD_NLOAD_MAX);
1159         }
1160     }
1161 }
1162 #endif
1163
1164 void dd_setup_dlb_resource_sharing(t_commrec            *cr,
1165                                    int                   gpu_id)
1166 {
1167 #if GMX_MPI
1168     int           physicalnode_id_hash;
1169     gmx_domdec_t *dd;
1170     MPI_Comm      mpi_comm_pp_physicalnode;
1171
1172     if (!thisRankHasDuty(cr, DUTY_PP) || gpu_id < 0)
1173     {
1174         /* Only ranks with short-ranged tasks (currently) use GPUs.
1175          * If we don't have GPUs assigned, there are no resources to share.
1176          */
1177         return;
1178     }
1179
1180     physicalnode_id_hash = gmx_physicalnode_id_hash();
1181
1182     dd = cr->dd;
1183
1184     if (debug)
1185     {
1186         fprintf(debug, "dd_setup_dd_dlb_gpu_sharing:\n");
1187         fprintf(debug, "DD PP rank %d physical node hash %d gpu_id %d\n",
1188                 dd->rank, physicalnode_id_hash, gpu_id);
1189     }
1190     /* Split the PP communicator over the physical nodes */
1191     /* TODO: See if we should store this (before), as it's also used for
1192      * for the nodecomm summation.
1193      */
1194     // TODO PhysicalNodeCommunicator could be extended/used to handle
1195     // the need for per-node per-group communicators.
1196     MPI_Comm_split(dd->mpi_comm_all, physicalnode_id_hash, dd->rank,
1197                    &mpi_comm_pp_physicalnode);
1198     MPI_Comm_split(mpi_comm_pp_physicalnode, gpu_id, dd->rank,
1199                    &dd->comm->mpi_comm_gpu_shared);
1200     MPI_Comm_free(&mpi_comm_pp_physicalnode);
1201     MPI_Comm_size(dd->comm->mpi_comm_gpu_shared, &dd->comm->nrank_gpu_shared);
1202
1203     if (debug)
1204     {
1205         fprintf(debug, "nrank_gpu_shared %d\n", dd->comm->nrank_gpu_shared);
1206     }
1207
1208     /* Note that some ranks could share a GPU, while others don't */
1209
1210     if (dd->comm->nrank_gpu_shared == 1)
1211     {
1212         MPI_Comm_free(&dd->comm->mpi_comm_gpu_shared);
1213     }
1214 #else
1215     GMX_UNUSED_VALUE(cr);
1216     GMX_UNUSED_VALUE(gpu_id);
1217 #endif
1218 }
1219
1220 static void make_load_communicators(gmx_domdec_t gmx_unused *dd)
1221 {
1222 #if GMX_MPI
1223     int  dim0, dim1, i, j;
1224     ivec loc;
1225
1226     if (debug)
1227     {
1228         fprintf(debug, "Making load communicators\n");
1229     }
1230
1231     dd->comm->load = new domdec_load_t[std::max(dd->ndim, 1)];
1232     snew(dd->comm->mpi_comm_load, std::max(dd->ndim, 1));
1233
1234     if (dd->ndim == 0)
1235     {
1236         return;
1237     }
1238
1239     clear_ivec(loc);
1240     make_load_communicator(dd, 0, loc);
1241     if (dd->ndim > 1)
1242     {
1243         dim0 = dd->dim[0];
1244         for (i = 0; i < dd->nc[dim0]; i++)
1245         {
1246             loc[dim0] = i;
1247             make_load_communicator(dd, 1, loc);
1248         }
1249     }
1250     if (dd->ndim > 2)
1251     {
1252         dim0 = dd->dim[0];
1253         for (i = 0; i < dd->nc[dim0]; i++)
1254         {
1255             loc[dim0] = i;
1256             dim1      = dd->dim[1];
1257             for (j = 0; j < dd->nc[dim1]; j++)
1258             {
1259                 loc[dim1] = j;
1260                 make_load_communicator(dd, 2, loc);
1261             }
1262         }
1263     }
1264
1265     if (debug)
1266     {
1267         fprintf(debug, "Finished making load communicators\n");
1268     }
1269 #endif
1270 }
1271
1272 /*! \brief Sets up the relation between neighboring domains and zones */
1273 static void setup_neighbor_relations(gmx_domdec_t *dd)
1274 {
1275     int                     d, dim, i, j, m;
1276     ivec                    tmp, s;
1277     gmx_domdec_zones_t     *zones;
1278     gmx_domdec_ns_ranges_t *izone;
1279     GMX_ASSERT(dd->ndim >= 0, "Must have non-negative number of dimensions for DD");
1280
1281     for (d = 0; d < dd->ndim; d++)
1282     {
1283         dim = dd->dim[d];
1284         copy_ivec(dd->ci, tmp);
1285         tmp[dim]           = (tmp[dim] + 1) % dd->nc[dim];
1286         dd->neighbor[d][0] = ddcoord2ddnodeid(dd, tmp);
1287         copy_ivec(dd->ci, tmp);
1288         tmp[dim]           = (tmp[dim] - 1 + dd->nc[dim]) % dd->nc[dim];
1289         dd->neighbor[d][1] = ddcoord2ddnodeid(dd, tmp);
1290         if (debug)
1291         {
1292             fprintf(debug, "DD rank %d neighbor ranks in dir %d are + %d - %d\n",
1293                     dd->rank, dim,
1294                     dd->neighbor[d][0],
1295                     dd->neighbor[d][1]);
1296         }
1297     }
1298
1299     int nzone  = (1 << dd->ndim);
1300     GMX_ASSERT(dd->ndim < DIM, "Invalid number of dimensions");
1301     int nizone = (1 << std::max(dd->ndim - 1, 0));
1302     assert(nizone >= 1 && nizone <= DD_MAXIZONE);
1303
1304     zones = &dd->comm->zones;
1305
1306     for (i = 0; i < nzone; i++)
1307     {
1308         m = 0;
1309         clear_ivec(zones->shift[i]);
1310         for (d = 0; d < dd->ndim; d++)
1311         {
1312             zones->shift[i][dd->dim[d]] = dd_zo[i][m++];
1313         }
1314     }
1315
1316     zones->n = nzone;
1317     for (i = 0; i < nzone; i++)
1318     {
1319         for (d = 0; d < DIM; d++)
1320         {
1321             s[d] = dd->ci[d] - zones->shift[i][d];
1322             if (s[d] < 0)
1323             {
1324                 s[d] += dd->nc[d];
1325             }
1326             else if (s[d] >= dd->nc[d])
1327             {
1328                 s[d] -= dd->nc[d];
1329             }
1330         }
1331     }
1332     zones->nizone = nizone;
1333     for (i = 0; i < zones->nizone; i++)
1334     {
1335         assert(ddNonbondedZonePairRanges[i][0] == i);
1336
1337         izone     = &zones->izone[i];
1338         /* dd_zp3 is for 3D decomposition, for fewer dimensions use only
1339          * j-zones up to nzone.
1340          */
1341         izone->j0 = std::min(ddNonbondedZonePairRanges[i][1], nzone);
1342         izone->j1 = std::min(ddNonbondedZonePairRanges[i][2], nzone);
1343         for (dim = 0; dim < DIM; dim++)
1344         {
1345             if (dd->nc[dim] == 1)
1346             {
1347                 /* All shifts should be allowed */
1348                 izone->shift0[dim] = -1;
1349                 izone->shift1[dim] = 1;
1350             }
1351             else
1352             {
1353                 /* Determine the min/max j-zone shift wrt the i-zone */
1354                 izone->shift0[dim] = 1;
1355                 izone->shift1[dim] = -1;
1356                 for (j = izone->j0; j < izone->j1; j++)
1357                 {
1358                     int shift_diff = zones->shift[j][dim] - zones->shift[i][dim];
1359                     if (shift_diff < izone->shift0[dim])
1360                     {
1361                         izone->shift0[dim] = shift_diff;
1362                     }
1363                     if (shift_diff > izone->shift1[dim])
1364                     {
1365                         izone->shift1[dim] = shift_diff;
1366                     }
1367                 }
1368             }
1369         }
1370     }
1371
1372     if (!isDlbDisabled(dd->comm))
1373     {
1374         dd->comm->cellsizesWithDlb.resize(dd->ndim);
1375     }
1376
1377     if (dd->comm->ddSettings.recordLoad)
1378     {
1379         make_load_communicators(dd);
1380     }
1381 }
1382
1383 static void make_pp_communicator(const gmx::MDLogger  &mdlog,
1384                                  gmx_domdec_t         *dd,
1385                                  t_commrec gmx_unused *cr,
1386                                  bool gmx_unused       reorder)
1387 {
1388 #if GMX_MPI
1389     gmx_domdec_comm_t *comm        = dd->comm;
1390     DDRankSetup       &ddRankSetup = comm->ddRankSetup;
1391
1392     if (ddRankSetup.bCartesianPP)
1393     {
1394         /* Set up cartesian communication for the particle-particle part */
1395         GMX_LOG(mdlog.info).appendTextFormatted(
1396                 "Will use a Cartesian communicator: %d x %d x %d",
1397                 dd->nc[XX], dd->nc[YY], dd->nc[ZZ]);
1398
1399         ivec periods;
1400         for (int i = 0; i < DIM; i++)
1401         {
1402             periods[i] = TRUE;
1403         }
1404         MPI_Comm comm_cart;
1405         MPI_Cart_create(cr->mpi_comm_mygroup, DIM, dd->nc, periods, static_cast<int>(reorder),
1406                         &comm_cart);
1407         /* We overwrite the old communicator with the new cartesian one */
1408         cr->mpi_comm_mygroup = comm_cart;
1409     }
1410
1411     dd->mpi_comm_all = cr->mpi_comm_mygroup;
1412     MPI_Comm_rank(dd->mpi_comm_all, &dd->rank);
1413
1414     if (ddRankSetup.bCartesianPP_PME)
1415     {
1416         /* Since we want to use the original cartesian setup for sim,
1417          * and not the one after split, we need to make an index.
1418          */
1419         snew(ddRankSetup.ddindex2ddnodeid, dd->nnodes);
1420         ddRankSetup.ddindex2ddnodeid[dd_index(dd->nc, dd->ci)] = dd->rank;
1421         gmx_sumi(dd->nnodes, ddRankSetup.ddindex2ddnodeid, cr);
1422         /* Get the rank of the DD master,
1423          * above we made sure that the master node is a PP node.
1424          */
1425         int rank;
1426         if (MASTER(cr))
1427         {
1428             rank = dd->rank;
1429         }
1430         else
1431         {
1432             rank = 0;
1433         }
1434         MPI_Allreduce(&rank, &dd->masterrank, 1, MPI_INT, MPI_SUM, dd->mpi_comm_all);
1435     }
1436     else if (ddRankSetup.bCartesianPP)
1437     {
1438         if (cr->npmenodes == 0)
1439         {
1440             /* The PP communicator is also
1441              * the communicator for this simulation
1442              */
1443             cr->mpi_comm_mysim = cr->mpi_comm_mygroup;
1444         }
1445         cr->nodeid = dd->rank;
1446
1447         MPI_Cart_coords(dd->mpi_comm_all, dd->rank, DIM, dd->ci);
1448
1449         /* We need to make an index to go from the coordinates
1450          * to the nodeid of this simulation.
1451          */
1452         snew(ddRankSetup.ddindex2simnodeid, dd->nnodes);
1453         std::vector<int> buf(dd->nnodes);
1454         if (thisRankHasDuty(cr, DUTY_PP))
1455         {
1456             buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
1457         }
1458         /* Communicate the ddindex to simulation nodeid index */
1459         MPI_Allreduce(buf.data(), ddRankSetup.ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
1460                       cr->mpi_comm_mysim);
1461
1462         /* Determine the master coordinates and rank.
1463          * The DD master should be the same node as the master of this sim.
1464          */
1465         for (int i = 0; i < dd->nnodes; i++)
1466         {
1467             if (ddRankSetup.ddindex2simnodeid[i] == 0)
1468             {
1469                 ddindex2xyz(dd->nc, i, dd->master_ci);
1470                 MPI_Cart_rank(dd->mpi_comm_all, dd->master_ci, &dd->masterrank);
1471             }
1472         }
1473         if (debug)
1474         {
1475             fprintf(debug, "The master rank is %d\n", dd->masterrank);
1476         }
1477     }
1478     else
1479     {
1480         /* No Cartesian communicators */
1481         /* We use the rank in dd->comm->all as DD index */
1482         ddindex2xyz(dd->nc, dd->rank, dd->ci);
1483         /* The simulation master nodeid is 0, so the DD master rank is also 0 */
1484         dd->masterrank = 0;
1485         clear_ivec(dd->master_ci);
1486     }
1487 #endif
1488
1489     GMX_LOG(mdlog.info).appendTextFormatted(
1490             "Domain decomposition rank %d, coordinates %d %d %d\n",
1491             dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
1492     if (debug)
1493     {
1494         fprintf(debug,
1495                 "Domain decomposition rank %d, coordinates %d %d %d\n\n",
1496                 dd->rank, dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
1497     }
1498 }
1499
1500 static void receive_ddindex2simnodeid(gmx_domdec_t         *dd,
1501                                       t_commrec            *cr)
1502 {
1503 #if GMX_MPI
1504     DDRankSetup &ddRankSetup = dd->comm->ddRankSetup;
1505
1506     if (!ddRankSetup.bCartesianPP_PME && ddRankSetup.bCartesianPP)
1507     {
1508         int *buf;
1509         snew(ddRankSetup.ddindex2simnodeid, dd->nnodes);
1510         snew(buf, dd->nnodes);
1511         if (thisRankHasDuty(cr, DUTY_PP))
1512         {
1513             buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
1514         }
1515         /* Communicate the ddindex to simulation nodeid index */
1516         MPI_Allreduce(buf, ddRankSetup.ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
1517                       cr->mpi_comm_mysim);
1518         sfree(buf);
1519     }
1520 #else
1521     GMX_UNUSED_VALUE(dd);
1522     GMX_UNUSED_VALUE(cr);
1523 #endif
1524 }
1525
1526 static void split_communicator(const gmx::MDLogger    &mdlog,
1527                                t_commrec              *cr,
1528                                DdRankOrder gmx_unused  rankOrder,
1529                                bool gmx_unused         reorder,
1530                                DDRankSetup            *ddRankSetup,
1531                                ivec                    ddCellIndex)
1532 {
1533     const ivec &numDDCells    = ddRankSetup->numPPCells;
1534
1535     if (ddRankSetup->bCartesianPP)
1536     {
1537         const int numDDCellsTot = ddRankSetup->numPPRanks;
1538         bool      bDiv[DIM];
1539         for (int i = 1; i < DIM; i++)
1540         {
1541             bDiv[i] = ((cr->npmenodes*numDDCells[i]) % numDDCellsTot == 0);
1542         }
1543         if (bDiv[YY] || bDiv[ZZ])
1544         {
1545             ddRankSetup->bCartesianPP_PME = TRUE;
1546             /* If we have 2D PME decomposition, which is always in x+y,
1547              * we stack the PME only nodes in z.
1548              * Otherwise we choose the direction that provides the thinnest slab
1549              * of PME only nodes as this will have the least effect
1550              * on the PP communication.
1551              * But for the PME communication the opposite might be better.
1552              */
1553             if (bDiv[ZZ] && (ddRankSetup->npmenodes_y > 1 ||
1554                              !bDiv[YY] ||
1555                              numDDCells[YY] > numDDCells[ZZ]))
1556             {
1557                 ddRankSetup->cartpmedim = ZZ;
1558             }
1559             else
1560             {
1561                 ddRankSetup->cartpmedim = YY;
1562             }
1563             ddRankSetup->ntot[ddRankSetup->cartpmedim]
1564                 += (cr->npmenodes*numDDCells[ddRankSetup->cartpmedim])/numDDCellsTot;
1565         }
1566         else
1567         {
1568             GMX_LOG(mdlog.info).appendTextFormatted(
1569                     "Number of PME-only ranks (%d) is not a multiple of nx*ny (%d*%d) or nx*nz (%d*%d)",
1570                     cr->npmenodes, numDDCells[XX], numDDCells[YY], numDDCells[XX], numDDCells[ZZ]);
1571             GMX_LOG(mdlog.info).appendText("Will not use a Cartesian communicator for PP <-> PME\n");
1572         }
1573     }
1574
1575     if (ddRankSetup->bCartesianPP_PME)
1576     {
1577 #if GMX_MPI
1578         int  rank;
1579         ivec periods;
1580
1581         GMX_LOG(mdlog.info).appendTextFormatted(
1582                 "Will use a Cartesian communicator for PP <-> PME: %d x %d x %d",
1583                 ddRankSetup->ntot[XX], ddRankSetup->ntot[YY], ddRankSetup->ntot[ZZ]);
1584
1585         for (int i = 0; i < DIM; i++)
1586         {
1587             periods[i] = TRUE;
1588         }
1589         MPI_Comm comm_cart;
1590         MPI_Cart_create(cr->mpi_comm_mysim, DIM, ddRankSetup->ntot, periods, static_cast<int>(reorder),
1591                         &comm_cart);
1592         MPI_Comm_rank(comm_cart, &rank);
1593         if (MASTER(cr) && rank != 0)
1594         {
1595             gmx_fatal(FARGS, "MPI rank 0 was renumbered by MPI_Cart_create, we do not allow this");
1596         }
1597
1598         /* With this assigment we loose the link to the original communicator
1599          * which will usually be MPI_COMM_WORLD, unless have multisim.
1600          */
1601         cr->mpi_comm_mysim = comm_cart;
1602         cr->sim_nodeid     = rank;
1603
1604         MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, ddCellIndex);
1605
1606         GMX_LOG(mdlog.info).appendTextFormatted(
1607                 "Cartesian rank %d, coordinates %d %d %d\n",
1608                 cr->sim_nodeid, ddCellIndex[XX], ddCellIndex[YY], ddCellIndex[ZZ]);
1609
1610         if (ddCellIndex[ddRankSetup->cartpmedim] < numDDCells[ddRankSetup->cartpmedim])
1611         {
1612             cr->duty = DUTY_PP;
1613         }
1614         if (cr->npmenodes == 0 ||
1615             ddCellIndex[ddRankSetup->cartpmedim] >= numDDCells[ddRankSetup->cartpmedim])
1616         {
1617             cr->duty = DUTY_PME;
1618         }
1619
1620         /* Split the sim communicator into PP and PME only nodes */
1621         MPI_Comm_split(cr->mpi_comm_mysim,
1622                        getThisRankDuties(cr),
1623                        dd_index(ddRankSetup->ntot, ddCellIndex),
1624                        &cr->mpi_comm_mygroup);
1625 #endif
1626     }
1627     else
1628     {
1629         switch (rankOrder)
1630         {
1631             case DdRankOrder::pp_pme:
1632                 GMX_LOG(mdlog.info).appendText("Order of the ranks: PP first, PME last");
1633                 break;
1634             case DdRankOrder::interleave:
1635                 /* Interleave the PP-only and PME-only ranks */
1636                 GMX_LOG(mdlog.info).appendText("Interleaving PP and PME ranks");
1637                 ddRankSetup->pmenodes = dd_interleaved_pme_ranks(*ddRankSetup);
1638                 break;
1639             case DdRankOrder::cartesian:
1640                 break;
1641             default:
1642                 gmx_fatal(FARGS, "Invalid ddRankOrder=%d", static_cast<int>(rankOrder));
1643         }
1644
1645         if (dd_simnode2pmenode(*ddRankSetup, cr, cr->sim_nodeid) == -1)
1646         {
1647             cr->duty = DUTY_PME;
1648         }
1649         else
1650         {
1651             cr->duty = DUTY_PP;
1652         }
1653 #if GMX_MPI
1654         /* Split the sim communicator into PP and PME only nodes */
1655         MPI_Comm_split(cr->mpi_comm_mysim,
1656                        getThisRankDuties(cr),
1657                        cr->nodeid,
1658                        &cr->mpi_comm_mygroup);
1659         MPI_Comm_rank(cr->mpi_comm_mygroup, &cr->nodeid);
1660 #endif
1661     }
1662
1663     GMX_LOG(mdlog.info).appendTextFormatted(
1664             "This rank does only %s work.\n",
1665             thisRankHasDuty(cr, DUTY_PP) ? "particle-particle" : "PME-mesh");
1666 }
1667
1668 /*! \brief Makes the PP communicator and the PME communicator, when needed
1669  *
1670  * Updates \p ddRankSetup for Cartesian communication aspects.
1671  * For PP ranks, sets the DD PP cell index in \p ddCellIndex.
1672  */
1673 static void makeGroupCommunicators(const gmx::MDLogger &mdlog,
1674                                    const DDSettings    &ddSettings,
1675                                    const DDSetup       &ddSetup,
1676                                    const DdRankOrder    ddRankOrder,
1677                                    DDRankSetup         *ddRankSetup,
1678                                    t_commrec           *cr,
1679                                    ivec                 ddCellIndex)
1680 {
1681     /* Initially we set ntot to the number of PP cells,
1682      * This will be increased with PME cells when using Cartesian communicators.
1683      */
1684     copy_ivec(ddSetup.numDomains, ddRankSetup->ntot);
1685
1686     ddRankSetup->bCartesianPP     = (ddRankOrder == DdRankOrder::cartesian);
1687     ddRankSetup->bCartesianPP_PME = FALSE;
1688
1689     if (ddSetup.numPmeRanks > 0)
1690     {
1691         /* Split the communicator into a PP and PME part */
1692         split_communicator(mdlog, cr, ddRankOrder, ddSettings.useCartesianReorder,
1693                            ddRankSetup, ddCellIndex);
1694     }
1695     else
1696     {
1697         /* All nodes do PP and PME */
1698         /* We do not require separate communicators */
1699         cr->mpi_comm_mygroup = cr->mpi_comm_mysim;
1700     }
1701 }
1702
1703 /*! \brief For PP ranks, sets or makes the communicator
1704  *
1705  * For PME ranks get the rank id.
1706  * For PP only ranks, sets the PME-only rank.
1707  */
1708 static void setupGroupCommunication(const gmx::MDLogger &mdlog,
1709                                     const DDSettings    &ddSettings,
1710                                     t_commrec           *cr,
1711                                     gmx_domdec_t        *dd)
1712 {
1713     DDRankSetup &ddRankSetup = dd->comm->ddRankSetup;
1714
1715     if (thisRankHasDuty(cr, DUTY_PP))
1716     {
1717         /* Copy or make a new PP communicator */
1718
1719         /* We (possibly) reordered the nodes in split_communicator,
1720          * so it is no longer required in make_pp_communicator.
1721          */
1722         const bool useCartesianReorder =
1723             (ddSettings.useCartesianReorder &&
1724              !ddRankSetup.bCartesianPP_PME);
1725
1726         make_pp_communicator(mdlog, dd, cr, useCartesianReorder);
1727     }
1728     else
1729     {
1730         receive_ddindex2simnodeid(dd, cr);
1731     }
1732
1733     if (!thisRankHasDuty(cr, DUTY_PME))
1734     {
1735         /* Set up the commnuication to our PME node */
1736         dd->pme_nodeid           = dd_simnode2pmenode(ddRankSetup, cr, cr->sim_nodeid);
1737         dd->pme_receive_vir_ener = receive_vir_ener(dd, cr);
1738         if (debug)
1739         {
1740             fprintf(debug, "My pme_nodeid %d receive ener %s\n",
1741                     dd->pme_nodeid, gmx::boolToString(dd->pme_receive_vir_ener));
1742         }
1743     }
1744     else
1745     {
1746         dd->pme_nodeid = -1;
1747     }
1748
1749     /* We can not use DDMASTER(dd), because dd->masterrank is set later */
1750     if (MASTER(cr))
1751     {
1752         dd->ma = std::make_unique<AtomDistribution>(dd->nc,
1753                                                     dd->comm->cgs_gl.nr,
1754                                                     dd->comm->cgs_gl.index[dd->comm->cgs_gl.nr]);
1755     }
1756 }
1757
1758 static real *get_slb_frac(const gmx::MDLogger &mdlog,
1759                           const char *dir, int nc, const char *size_string)
1760 {
1761     real  *slb_frac, tot;
1762     int    i, n;
1763     double dbl;
1764
1765     slb_frac = nullptr;
1766     if (nc > 1 && size_string != nullptr)
1767     {
1768         GMX_LOG(mdlog.info).appendTextFormatted(
1769                 "Using static load balancing for the %s direction", dir);
1770         snew(slb_frac, nc);
1771         tot = 0;
1772         for (i = 0; i < nc; i++)
1773         {
1774             dbl = 0;
1775             sscanf(size_string, "%20lf%n", &dbl, &n);
1776             if (dbl == 0)
1777             {
1778                 gmx_fatal(FARGS, "Incorrect or not enough DD cell size entries for direction %s: '%s'", dir, size_string);
1779             }
1780             slb_frac[i]  = dbl;
1781             size_string += n;
1782             tot         += slb_frac[i];
1783         }
1784         /* Normalize */
1785         std::string relativeCellSizes = "Relative cell sizes:";
1786         for (i = 0; i < nc; i++)
1787         {
1788             slb_frac[i]       /= tot;
1789             relativeCellSizes += gmx::formatString(" %5.3f", slb_frac[i]);
1790         }
1791         GMX_LOG(mdlog.info).appendText(relativeCellSizes);
1792     }
1793
1794     return slb_frac;
1795 }
1796
1797 static int multi_body_bondeds_count(const gmx_mtop_t *mtop)
1798 {
1799     int                  n     = 0;
1800     gmx_mtop_ilistloop_t iloop = gmx_mtop_ilistloop_init(mtop);
1801     int                  nmol;
1802     while (const InteractionLists *ilists = gmx_mtop_ilistloop_next(iloop, &nmol))
1803     {
1804         for (auto &ilist : extractILists(*ilists, IF_BOND))
1805         {
1806             if (NRAL(ilist.functionType) >  2)
1807             {
1808                 n += nmol*(ilist.iatoms.size()/ilistStride(ilist));
1809             }
1810         }
1811     }
1812
1813     return n;
1814 }
1815
1816 static int dd_getenv(const gmx::MDLogger &mdlog,
1817                      const char *env_var, int def)
1818 {
1819     char *val;
1820     int   nst;
1821
1822     nst = def;
1823     val = getenv(env_var);
1824     if (val)
1825     {
1826         if (sscanf(val, "%20d", &nst) <= 0)
1827         {
1828             nst = 1;
1829         }
1830         GMX_LOG(mdlog.info).appendTextFormatted(
1831                 "Found env.var. %s = %s, using value %d",
1832                 env_var, val, nst);
1833     }
1834
1835     return nst;
1836 }
1837
1838 static void check_dd_restrictions(const gmx_domdec_t  *dd,
1839                                   const t_inputrec    *ir,
1840                                   const gmx::MDLogger &mdlog)
1841 {
1842     if (ir->ePBC == epbcSCREW &&
1843         (dd->nc[XX] == 1 || dd->nc[YY] > 1 || dd->nc[ZZ] > 1))
1844     {
1845         gmx_fatal(FARGS, "With pbc=%s can only do domain decomposition in the x-direction", epbc_names[ir->ePBC]);
1846     }
1847
1848     if (ir->ns_type == ensSIMPLE)
1849     {
1850         gmx_fatal(FARGS, "Domain decomposition does not support simple neighbor searching, use grid searching or run with one MPI rank");
1851     }
1852
1853     if (ir->nstlist == 0)
1854     {
1855         gmx_fatal(FARGS, "Domain decomposition does not work with nstlist=0");
1856     }
1857
1858     if (ir->comm_mode == ecmANGULAR && ir->ePBC != epbcNONE)
1859     {
1860         GMX_LOG(mdlog.warning).appendText("comm-mode angular will give incorrect results when the comm group partially crosses a periodic boundary");
1861     }
1862 }
1863
1864 static real average_cellsize_min(const gmx_ddbox_t &ddbox,
1865                                  const ivec         numDomains)
1866 {
1867     real r = ddbox.box_size[XX];
1868     for (int d = 0; d < DIM; d++)
1869     {
1870         if (numDomains[d] > 1)
1871         {
1872             /* Check using the initial average cell size */
1873             r = std::min(r, ddbox.box_size[d]*ddbox.skew_fac[d]/numDomains[d]);
1874         }
1875     }
1876
1877     return r;
1878 }
1879
1880 /*! \brief Depending on the DLB initial value return the DLB switched off state or issue an error.
1881  */
1882 static DlbState forceDlbOffOrBail(DlbState             cmdlineDlbState,
1883                                   const std::string   &reasonStr,
1884                                   const gmx::MDLogger &mdlog)
1885 {
1886     std::string dlbNotSupportedErr  = "Dynamic load balancing requested, but ";
1887     std::string dlbDisableNote      = "NOTE: disabling dynamic load balancing as ";
1888
1889     if (cmdlineDlbState == DlbState::onUser)
1890     {
1891         gmx_fatal(FARGS, "%s", (dlbNotSupportedErr + reasonStr).c_str());
1892     }
1893     else if (cmdlineDlbState == DlbState::offCanTurnOn)
1894     {
1895         GMX_LOG(mdlog.info).appendText(dlbDisableNote + reasonStr);
1896     }
1897     return DlbState::offForever;
1898 }
1899
1900 /*! \brief Return the dynamic load balancer's initial state based on initial conditions and user inputs.
1901  *
1902  * This function parses the parameters of "-dlb" command line option setting
1903  * corresponding state values. Then it checks the consistency of the determined
1904  * state with other run parameters and settings. As a result, the initial state
1905  * may be altered or an error may be thrown if incompatibility of options is detected.
1906  *
1907  * \param [in] mdlog       Logger.
1908  * \param [in] dlbOption   Enum value for the DLB option.
1909  * \param [in] bRecordLoad True if the load balancer is recording load information.
1910  * \param [in] mdrunOptions  Options for mdrun.
1911  * \param [in] ir          Pointer mdrun to input parameters.
1912  * \returns                DLB initial/startup state.
1913  */
1914 static DlbState determineInitialDlbState(const gmx::MDLogger &mdlog,
1915                                          DlbOption dlbOption, gmx_bool bRecordLoad,
1916                                          const gmx::MdrunOptions &mdrunOptions,
1917                                          const t_inputrec *ir)
1918 {
1919     DlbState dlbState = DlbState::offCanTurnOn;
1920
1921     switch (dlbOption)
1922     {
1923         case DlbOption::turnOnWhenUseful: dlbState = DlbState::offCanTurnOn; break;
1924         case DlbOption::no:               dlbState = DlbState::offUser;      break;
1925         case DlbOption::yes:              dlbState = DlbState::onUser;       break;
1926         default: gmx_incons("Invalid dlbOption enum value");
1927     }
1928
1929     /* Reruns don't support DLB: bail or override auto mode */
1930     if (mdrunOptions.rerun)
1931     {
1932         std::string reasonStr = "it is not supported in reruns.";
1933         return forceDlbOffOrBail(dlbState, reasonStr, mdlog);
1934     }
1935
1936     /* Unsupported integrators */
1937     if (!EI_DYNAMICS(ir->eI))
1938     {
1939         auto reasonStr = gmx::formatString("it is only supported with dynamics, not with integrator '%s'.", EI(ir->eI));
1940         return forceDlbOffOrBail(dlbState, reasonStr, mdlog);
1941     }
1942
1943     /* Without cycle counters we can't time work to balance on */
1944     if (!bRecordLoad)
1945     {
1946         std::string reasonStr = "cycle counters unsupported or not enabled in the operating system kernel.";
1947         return forceDlbOffOrBail(dlbState, reasonStr, mdlog);
1948     }
1949
1950     if (mdrunOptions.reproducible)
1951     {
1952         std::string reasonStr = "you started a reproducible run.";
1953         switch (dlbState)
1954         {
1955             case DlbState::offUser:
1956                 break;
1957             case DlbState::offForever:
1958                 GMX_RELEASE_ASSERT(false, "DlbState::offForever is not a valid initial state");
1959                 break;
1960             case DlbState::offCanTurnOn:
1961                 return forceDlbOffOrBail(dlbState, reasonStr, mdlog);
1962             case DlbState::onCanTurnOff:
1963                 GMX_RELEASE_ASSERT(false, "DlbState::offCanTurnOff is not a valid initial state");
1964                 break;
1965             case DlbState::onUser:
1966                 return forceDlbOffOrBail(dlbState, reasonStr + " In load balanced runs binary reproducibility cannot be ensured.", mdlog);
1967             default:
1968                 gmx_fatal(FARGS, "Death horror: undefined case (%d) for load balancing choice", static_cast<int>(dlbState));
1969         }
1970     }
1971
1972     return dlbState;
1973 }
1974
1975 /* Sets the order of the DD dimensions, returns the number of DD dimensions */
1976 static int set_dd_dim(const ivec        numDDCells,
1977                       const DDSettings &ddSettings,
1978                       ivec              dims)
1979 {
1980     int ndim = 0;
1981     if (ddSettings.useDDOrderZYX)
1982     {
1983         /* Decomposition order z,y,x */
1984         for (int dim = DIM - 1; dim >= 0; dim--)
1985         {
1986             if (numDDCells[dim] > 1)
1987             {
1988                 dims[ndim++] = dim;
1989             }
1990         }
1991     }
1992     else
1993     {
1994         /* Decomposition order x,y,z */
1995         for (int dim = 0; dim < DIM; dim++)
1996         {
1997             if (numDDCells[dim] > 1)
1998             {
1999                 dims[ndim++] = dim;
2000             }
2001         }
2002     }
2003
2004     if (ndim == 0)
2005     {
2006         /* Set dim[0] to avoid extra checks on ndim in several places */
2007         dims[0] = XX;
2008     }
2009
2010     return ndim;
2011 }
2012
2013 static gmx_domdec_comm_t *init_dd_comm()
2014 {
2015     gmx_domdec_comm_t *comm = new gmx_domdec_comm_t;
2016
2017     comm->n_load_have      = 0;
2018     comm->n_load_collect   = 0;
2019
2020     comm->haveTurnedOffDlb = false;
2021
2022     for (int i = 0; i < static_cast<int>(DDAtomRanges::Type::Number); i++)
2023     {
2024         comm->sum_nat[i] = 0;
2025     }
2026     comm->ndecomp   = 0;
2027     comm->nload     = 0;
2028     comm->load_step = 0;
2029     comm->load_sum  = 0;
2030     comm->load_max  = 0;
2031     clear_ivec(comm->load_lim);
2032     comm->load_mdf  = 0;
2033     comm->load_pme  = 0;
2034
2035     /* This should be replaced by a unique pointer */
2036     comm->balanceRegion = ddBalanceRegionAllocate();
2037
2038     return comm;
2039 }
2040
2041 /* Returns whether mtop contains constraints and/or vsites */
2042 static bool systemHasConstraintsOrVsites(const gmx_mtop_t &mtop)
2043 {
2044     auto ilistLoop = gmx_mtop_ilistloop_init(mtop);
2045     int  nmol;
2046     while (const InteractionLists *ilists = gmx_mtop_ilistloop_next(ilistLoop, &nmol))
2047     {
2048         if (!extractILists(*ilists, IF_CONSTRAINT | IF_VSITE).empty())
2049         {
2050             return true;
2051         }
2052     }
2053
2054     return false;
2055 }
2056
2057 static void setupUpdateGroups(const gmx::MDLogger &mdlog,
2058                               const gmx_mtop_t    &mtop,
2059                               const t_inputrec    &inputrec,
2060                               const real           cutoffMargin,
2061                               DDSystemInfo        *systemInfo)
2062 {
2063     /* When we have constraints and/or vsites, it is beneficial to use
2064      * update groups (when possible) to allow independent update of groups.
2065      */
2066     if (!systemHasConstraintsOrVsites(mtop))
2067     {
2068         /* No constraints or vsites, atoms can be updated independently */
2069         return;
2070     }
2071
2072     systemInfo->updateGroupingPerMoleculetype = gmx::makeUpdateGroups(mtop);
2073     systemInfo->useUpdateGroups               =
2074         (!systemInfo->updateGroupingPerMoleculetype.empty() &&
2075          getenv("GMX_NO_UPDATEGROUPS") == nullptr);
2076
2077     if (systemInfo->useUpdateGroups)
2078     {
2079         int numUpdateGroups = 0;
2080         for (const auto &molblock : mtop.molblock)
2081         {
2082             numUpdateGroups += molblock.nmol*systemInfo->updateGroupingPerMoleculetype[molblock.type].numBlocks();
2083         }
2084
2085         systemInfo->maxUpdateGroupRadius =
2086             computeMaxUpdateGroupRadius(mtop,
2087                                         systemInfo->updateGroupingPerMoleculetype,
2088                                         maxReferenceTemperature(inputrec));
2089
2090         /* To use update groups, the large domain-to-domain cutoff distance
2091          * should be compatible with the box size.
2092          */
2093         systemInfo->useUpdateGroups = (atomToAtomIntoDomainToDomainCutoff(*systemInfo, 0) < cutoffMargin);
2094
2095         if (systemInfo->useUpdateGroups)
2096         {
2097             GMX_LOG(mdlog.info).appendTextFormatted(
2098                     "Using update groups, nr %d, average size %.1f atoms, max. radius %.3f nm\n",
2099                     numUpdateGroups,
2100                     mtop.natoms/static_cast<double>(numUpdateGroups),
2101                     systemInfo->maxUpdateGroupRadius);
2102         }
2103         else
2104         {
2105             GMX_LOG(mdlog.info).appendTextFormatted("The combination of rlist and box size prohibits the use of update groups\n");
2106             systemInfo->updateGroupingPerMoleculetype.clear();
2107         }
2108     }
2109 }
2110
2111 UnitCellInfo::UnitCellInfo(const t_inputrec &ir) :
2112     npbcdim(ePBC2npbcdim(ir.ePBC)),
2113     numBoundedDimensions(inputrec2nboundeddim(&ir)),
2114     ddBoxIsDynamic(numBoundedDimensions < DIM || inputrecDynamicBox(&ir)),
2115     haveScrewPBC(ir.ePBC == epbcSCREW)
2116 {
2117 }
2118
2119 /*! \brief Generate the simulation system information */
2120 static DDSystemInfo
2121 getSystemInfo(const gmx::MDLogger           &mdlog,
2122               t_commrec                     *cr,
2123               const DomdecOptions           &options,
2124               const gmx_mtop_t              *mtop,
2125               const t_inputrec              *ir,
2126               const matrix                   box,
2127               gmx::ArrayRef<const gmx::RVec> xGlobal)
2128 {
2129     const real         tenPercentMargin = 1.1;
2130
2131     DDSystemInfo       systemInfo;
2132
2133     /* We need to decide on update groups early, as this affects communication distances */
2134     systemInfo.useUpdateGroups = false;
2135     if (ir->cutoff_scheme == ecutsVERLET)
2136     {
2137         real cutoffMargin = std::sqrt(max_cutoff2(ir->ePBC, box)) - ir->rlist;
2138         setupUpdateGroups(mdlog, *mtop, *ir, cutoffMargin, &systemInfo);
2139     }
2140
2141     // TODO: Check whether all bondeds are within update groups
2142     systemInfo.haveInterDomainBondeds          = (mtop->natoms > gmx_mtop_num_molecules(*mtop) ||
2143                                                   mtop->bIntermolecularInteractions);
2144     systemInfo.haveInterDomainMultiBodyBondeds = (multi_body_bondeds_count(mtop) > 0);
2145
2146     if (systemInfo.useUpdateGroups)
2147     {
2148         systemInfo.haveSplitConstraints = false;
2149         systemInfo.haveSplitSettles     = false;
2150     }
2151     else
2152     {
2153         systemInfo.haveSplitConstraints = gmx::inter_charge_group_constraints(*mtop);
2154         systemInfo.haveSplitSettles     = gmx::inter_charge_group_settles(*mtop);
2155     }
2156
2157     if (ir->rlist == 0)
2158     {
2159         /* Set the cut-off to some very large value,
2160          * so we don't need if statements everywhere in the code.
2161          * We use sqrt, since the cut-off is squared in some places.
2162          */
2163         systemInfo.cutoff = GMX_CUTOFF_INF;
2164     }
2165     else
2166     {
2167         systemInfo.cutoff = atomToAtomIntoDomainToDomainCutoff(systemInfo, ir->rlist);
2168     }
2169     systemInfo.minCutoffForMultiBody = 0;
2170
2171     /* Determine the minimum cell size limit, affected by many factors */
2172     systemInfo.cellsizeLimit             = 0;
2173     systemInfo.filterBondedCommunication = false;
2174
2175     /* We do not allow home atoms to move beyond the neighboring domain
2176      * between domain decomposition steps, which limits the cell size.
2177      * Get an estimate of cell size limit due to atom displacement.
2178      * In most cases this is a large overestimate, because it assumes
2179      * non-interaction atoms.
2180      * We set the chance to 1 in a trillion steps.
2181      */
2182     constexpr real c_chanceThatAtomMovesBeyondDomain = 1e-12;
2183     const real     limitForAtomDisplacement          =
2184         minCellSizeForAtomDisplacement(*mtop, *ir,
2185                                        systemInfo.updateGroupingPerMoleculetype,
2186                                        c_chanceThatAtomMovesBeyondDomain);
2187     GMX_LOG(mdlog.info).appendTextFormatted(
2188             "Minimum cell size due to atom displacement: %.3f nm",
2189             limitForAtomDisplacement);
2190
2191     systemInfo.cellsizeLimit = std::max(systemInfo.cellsizeLimit,
2192                                         limitForAtomDisplacement);
2193
2194     /* TODO: PME decomposition currently requires atoms not to be more than
2195      *       2/3 of comm->cutoff, which is >=rlist, outside of their domain.
2196      *       In nearly all cases, limitForAtomDisplacement will be smaller
2197      *       than 2/3*rlist, so the PME requirement is satisfied.
2198      *       But it would be better for both correctness and performance
2199      *       to use limitForAtomDisplacement instead of 2/3*comm->cutoff.
2200      *       Note that we would need to improve the pairlist buffer case.
2201      */
2202
2203     if (systemInfo.haveInterDomainBondeds)
2204     {
2205         if (options.minimumCommunicationRange > 0)
2206         {
2207             systemInfo.minCutoffForMultiBody =
2208                 atomToAtomIntoDomainToDomainCutoff(systemInfo, options.minimumCommunicationRange);
2209             if (options.useBondedCommunication)
2210             {
2211                 systemInfo.filterBondedCommunication = (systemInfo.minCutoffForMultiBody > systemInfo.cutoff);
2212             }
2213             else
2214             {
2215                 systemInfo.cutoff = std::max(systemInfo.cutoff,
2216                                              systemInfo.minCutoffForMultiBody);
2217             }
2218         }
2219         else if (ir->bPeriodicMols)
2220         {
2221             /* Can not easily determine the required cut-off */
2222             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.");
2223             systemInfo.minCutoffForMultiBody = systemInfo.cutoff/2;
2224         }
2225         else
2226         {
2227             real r_2b, r_mb;
2228
2229             if (MASTER(cr))
2230             {
2231                 dd_bonded_cg_distance(mdlog, mtop, ir, as_rvec_array(xGlobal.data()), box,
2232                                       options.checkBondedInteractions,
2233                                       &r_2b, &r_mb);
2234             }
2235             gmx_bcast(sizeof(r_2b), &r_2b, cr);
2236             gmx_bcast(sizeof(r_mb), &r_mb, cr);
2237
2238             /* We use an initial margin of 10% for the minimum cell size,
2239              * except when we are just below the non-bonded cut-off.
2240              */
2241             if (options.useBondedCommunication)
2242             {
2243                 if (std::max(r_2b, r_mb) > systemInfo.cutoff)
2244                 {
2245                     const real r_bonded              = std::max(r_2b, r_mb);
2246                     systemInfo.minCutoffForMultiBody = tenPercentMargin*r_bonded;
2247                     /* This is the (only) place where we turn on the filtering */
2248                     systemInfo.filterBondedCommunication = true;
2249                 }
2250                 else
2251                 {
2252                     const real r_bonded              = r_mb;
2253                     systemInfo.minCutoffForMultiBody = std::min(tenPercentMargin*r_bonded,
2254                                                                 systemInfo.cutoff);
2255                 }
2256                 /* We determine cutoff_mbody later */
2257                 systemInfo.increaseMultiBodyCutoff = true;
2258             }
2259             else
2260             {
2261                 /* No special bonded communication,
2262                  * simply increase the DD cut-off.
2263                  */
2264                 systemInfo.minCutoffForMultiBody = tenPercentMargin*std::max(r_2b, r_mb);
2265                 systemInfo.cutoff                = std::max(systemInfo.cutoff,
2266                                                             systemInfo.minCutoffForMultiBody);
2267             }
2268         }
2269         GMX_LOG(mdlog.info).appendTextFormatted(
2270                 "Minimum cell size due to bonded interactions: %.3f nm",
2271                 systemInfo.minCutoffForMultiBody);
2272
2273         systemInfo.cellsizeLimit = std::max(systemInfo.cellsizeLimit,
2274                                             systemInfo.minCutoffForMultiBody);
2275     }
2276
2277     systemInfo.constraintCommunicationRange = 0;
2278     if (systemInfo.haveSplitConstraints && options.constraintCommunicationRange <= 0)
2279     {
2280         /* There is a cell size limit due to the constraints (P-LINCS) */
2281         systemInfo.constraintCommunicationRange = gmx::constr_r_max(mdlog, mtop, ir);
2282         GMX_LOG(mdlog.info).appendTextFormatted(
2283                 "Estimated maximum distance required for P-LINCS: %.3f nm",
2284                 systemInfo.constraintCommunicationRange);
2285         if (systemInfo.constraintCommunicationRange > systemInfo.cellsizeLimit)
2286         {
2287             GMX_LOG(mdlog.info).appendText("This distance will limit the DD cell size, you can override this with -rcon");
2288         }
2289     }
2290     else if (options.constraintCommunicationRange > 0)
2291     {
2292         /* Here we do not check for dd->splitConstraints.
2293          * because one can also set a cell size limit for virtual sites only
2294          * and at this point we don't know yet if there are intercg v-sites.
2295          */
2296         GMX_LOG(mdlog.info).appendTextFormatted(
2297                 "User supplied maximum distance required for P-LINCS: %.3f nm",
2298                 options.constraintCommunicationRange);
2299         systemInfo.constraintCommunicationRange = options.constraintCommunicationRange;
2300     }
2301     systemInfo.cellsizeLimit = std::max(systemInfo.cellsizeLimit,
2302                                         systemInfo.constraintCommunicationRange);
2303
2304     return systemInfo;
2305 }
2306
2307 /*! \brief Set the cell size and interaction limits, as well as the DD grid
2308  *
2309  * Also computes the initial ddbox.
2310  */
2311 static DDSetup
2312 getDDSetup(const gmx::MDLogger           &mdlog,
2313            t_commrec                     *cr,
2314            const DomdecOptions           &options,
2315            const DDSettings              &ddSettings,
2316            const DDSystemInfo            &systemInfo,
2317            const gmx_mtop_t              *mtop,
2318            const t_inputrec              *ir,
2319            const matrix                   box,
2320            gmx::ArrayRef<const gmx::RVec> xGlobal,
2321            gmx_ddbox_t                   *ddbox)
2322 {
2323     DDSetup     ddSetup;
2324
2325     if (options.numCells[XX] > 0)
2326     {
2327         copy_ivec(options.numCells, ddSetup.numDomains);
2328         set_ddbox_cr(*cr, &ddSetup.numDomains, *ir, box, xGlobal, ddbox);
2329
2330         if (options.numPmeRanks >= 0)
2331         {
2332             ddSetup.numPmeRanks = options.numPmeRanks;
2333         }
2334         else
2335         {
2336             /* When the DD grid is set explicitly and -npme is set to auto,
2337              * don't use PME ranks. We check later if the DD grid is
2338              * compatible with the total number of ranks.
2339              */
2340             ddSetup.numPmeRanks = 0;
2341         }
2342     }
2343     else
2344     {
2345         set_ddbox_cr(*cr, nullptr, *ir, box, xGlobal, ddbox);
2346
2347         /* We need to choose the optimal DD grid and possibly PME nodes */
2348         ddSetup =
2349             dd_choose_grid(mdlog, cr, ir, mtop, box, ddbox,
2350                            options.numPmeRanks,
2351                            !isDlbDisabled(ddSettings.initialDlbState),
2352                            options.dlbScaling,
2353                            systemInfo);
2354
2355         if (ddSetup.numDomains[XX] == 0)
2356         {
2357             char     buf[STRLEN];
2358             gmx_bool bC = (systemInfo.haveSplitConstraints &&
2359                            systemInfo.constraintCommunicationRange > systemInfo.minCutoffForMultiBody);
2360             sprintf(buf, "Change the number of ranks or mdrun option %s%s%s",
2361                     !bC ? "-rdd" : "-rcon",
2362                     ddSettings.initialDlbState != DlbState::offUser ? " or -dds" : "",
2363                     bC ? " or your LINCS settings" : "");
2364
2365             gmx_fatal_collective(FARGS, cr->mpi_comm_mysim, MASTER(cr),
2366                                  "There is no domain decomposition for %d ranks that is compatible with the given box and a minimum cell size of %g nm\n"
2367                                  "%s\n"
2368                                  "Look in the log file for details on the domain decomposition",
2369                                  cr->nnodes - ddSetup.numPmeRanks,
2370                                  ddSetup.cellsizeLimit,
2371                                  buf);
2372         }
2373     }
2374
2375     const real acs = average_cellsize_min(*ddbox, ddSetup.numDomains);
2376     if (acs < systemInfo.cellsizeLimit)
2377     {
2378         if (options.numCells[XX] <= 0)
2379         {
2380             GMX_RELEASE_ASSERT(false, "dd_choose_grid() should return a grid that satisfies the cell size limits");
2381         }
2382         else
2383         {
2384             gmx_fatal_collective(FARGS, cr->mpi_comm_mysim, MASTER(cr),
2385                                  "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",
2386                                  acs, systemInfo.cellsizeLimit);
2387         }
2388     }
2389
2390     const int numPPRanks = ddSetup.numDomains[XX]*ddSetup.numDomains[YY]*ddSetup.numDomains[ZZ];
2391     if (cr->nnodes - numPPRanks != ddSetup.numPmeRanks)
2392     {
2393         gmx_fatal_collective(FARGS, cr->mpi_comm_mysim, MASTER(cr),
2394                              "The size of the domain decomposition grid (%d) does not match the number of ranks (%d). The total number of ranks is %d",
2395                              numPPRanks, cr->nnodes - ddSetup.numPmeRanks, cr->nnodes);
2396     }
2397     if (ddSetup.numPmeRanks > numPPRanks)
2398     {
2399         gmx_fatal_collective(FARGS, cr->mpi_comm_mysim, MASTER(cr),
2400                              "The number of separate PME ranks (%d) is larger than the number of PP ranks (%d), this is not supported.", ddSetup.numPmeRanks, numPPRanks);
2401     }
2402
2403     ddSetup.numDDDimensions = set_dd_dim(ddSetup.numDomains, ddSettings,
2404                                          ddSetup.ddDimensions);
2405
2406     return ddSetup;
2407 }
2408
2409 /*! \brief Set the cell size and interaction limits, as well as the DD grid */
2410 static DDRankSetup
2411 getDDRankSetup(const gmx::MDLogger &mdlog,
2412                t_commrec           *cr,
2413                const DDSetup       &ddSetup,
2414                const t_inputrec    &ir)
2415 {
2416     GMX_LOG(mdlog.info).appendTextFormatted(
2417             "Domain decomposition grid %d x %d x %d, separate PME ranks %d",
2418             ddSetup.numDomains[XX], ddSetup.numDomains[YY], ddSetup.numDomains[ZZ],
2419             ddSetup.numPmeRanks);
2420
2421     DDRankSetup ddRankSetup;
2422
2423     ddRankSetup.numPPRanks = cr->nnodes - cr->npmenodes;
2424     copy_ivec(ddSetup.numDomains, ddRankSetup.numPPCells);
2425
2426     if (cr->npmenodes > 0)
2427     {
2428         ddRankSetup.npmenodes = cr->npmenodes;
2429     }
2430     else
2431     {
2432         ddRankSetup.npmenodes = ddSetup.numDomains[XX]*ddSetup.numDomains[YY]*ddSetup.numDomains[ZZ];
2433     }
2434
2435     if (EEL_PME(ir.coulombtype) || EVDW_PME(ir.vdwtype))
2436     {
2437         /* The following choices should match those
2438          * in comm_cost_est in domdec_setup.c.
2439          * Note that here the checks have to take into account
2440          * that the decomposition might occur in a different order than xyz
2441          * (for instance through the env.var. GMX_DD_ORDER_ZYX),
2442          * in which case they will not match those in comm_cost_est,
2443          * but since that is mainly for testing purposes that's fine.
2444          */
2445         if (ddSetup.numDDDimensions >= 2 &&
2446             ddSetup.ddDimensions[0] == XX &&
2447             ddSetup.ddDimensions[1] == YY &&
2448             ddRankSetup.npmenodes > ddSetup.numDomains[XX] &&
2449             ddRankSetup.npmenodes % ddSetup.numDomains[XX] == 0 &&
2450             getenv("GMX_PMEONEDD") == nullptr)
2451         {
2452             ddRankSetup.npmedecompdim = 2;
2453             ddRankSetup.npmenodes_x   = ddSetup.numDomains[XX];
2454             ddRankSetup.npmenodes_y   = ddRankSetup.npmenodes/ddRankSetup.npmenodes_x;
2455         }
2456         else
2457         {
2458             /* In case nc is 1 in both x and y we could still choose to
2459              * decompose pme in y instead of x, but we use x for simplicity.
2460              */
2461             ddRankSetup.npmedecompdim = 1;
2462             if (ddSetup.ddDimensions[0] == YY)
2463             {
2464                 ddRankSetup.npmenodes_x = 1;
2465                 ddRankSetup.npmenodes_y = ddRankSetup.npmenodes;
2466             }
2467             else
2468             {
2469                 ddRankSetup.npmenodes_x = ddRankSetup.npmenodes;
2470                 ddRankSetup.npmenodes_y = 1;
2471             }
2472         }
2473         GMX_LOG(mdlog.info).appendTextFormatted(
2474                 "PME domain decomposition: %d x %d x %d",
2475                 ddRankSetup.npmenodes_x, ddRankSetup.npmenodes_y, 1);
2476     }
2477     else
2478     {
2479         ddRankSetup.npmedecompdim = 0;
2480         ddRankSetup.npmenodes_x   = 0;
2481         ddRankSetup.npmenodes_y   = 0;
2482     }
2483
2484     return ddRankSetup;
2485 }
2486
2487 /*! \brief Set the cell size and interaction limits */
2488 static void set_dd_limits(const gmx::MDLogger &mdlog,
2489                           t_commrec *cr, gmx_domdec_t *dd,
2490                           const DomdecOptions &options,
2491                           const DDSettings &ddSettings,
2492                           const DDSystemInfo &systemInfo,
2493                           const DDSetup &ddSetup,
2494                           const gmx_mtop_t *mtop,
2495                           const t_inputrec *ir,
2496                           const gmx_ddbox_t &ddbox)
2497 {
2498     gmx_domdec_comm_t *comm = dd->comm;
2499     comm->ddSettings        = ddSettings;
2500
2501     /* Initialize to GPU share count to 0, might change later */
2502     comm->nrank_gpu_shared = 0;
2503
2504     comm->dlbState         = comm->ddSettings.initialDlbState;
2505     dd_dlb_set_should_check_whether_to_turn_dlb_on(dd, TRUE);
2506     /* To consider turning DLB on after 2*nstlist steps we need to check
2507      * at partitioning count 3. Thus we need to increase the first count by 2.
2508      */
2509     comm->ddPartioningCountFirstDlbOff += 2;
2510
2511     comm->bPMELoadBalDLBLimits          = FALSE;
2512
2513     /* Allocate the charge group/atom sorting struct */
2514     comm->sort = std::make_unique<gmx_domdec_sort_t>();
2515
2516     comm->systemInfo = systemInfo;
2517
2518     if (systemInfo.useUpdateGroups)
2519     {
2520         /* Note: We would like to use dd->nnodes for the atom count estimate,
2521          *       but that is not yet available here. But this anyhow only
2522          *       affect performance up to the second dd_partition_system call.
2523          */
2524         const int homeAtomCountEstimate =  mtop->natoms/cr->nnodes;
2525         comm->updateGroupsCog =
2526             std::make_unique<gmx::UpdateGroupsCog>(*mtop,
2527                                                    systemInfo.updateGroupingPerMoleculetype,
2528                                                    maxReferenceTemperature(*ir),
2529                                                    homeAtomCountEstimate);
2530     }
2531
2532     comm->cgs_gl = gmx_mtop_global_cgs(mtop);
2533
2534     /* Set the DD setup given by ddSetup */
2535     cr->npmenodes = ddSetup.numPmeRanks;
2536     copy_ivec(ddSetup.numDomains, dd->nc);
2537     dd->ndim = ddSetup.numDDDimensions;
2538     copy_ivec(ddSetup.ddDimensions, dd->dim);
2539
2540     dd->nnodes = dd->nc[XX]*dd->nc[YY]*dd->nc[ZZ];
2541
2542     snew(comm->slb_frac, DIM);
2543     if (isDlbDisabled(comm))
2544     {
2545         comm->slb_frac[XX] = get_slb_frac(mdlog, "x", dd->nc[XX], options.cellSizeX);
2546         comm->slb_frac[YY] = get_slb_frac(mdlog, "y", dd->nc[YY], options.cellSizeY);
2547         comm->slb_frac[ZZ] = get_slb_frac(mdlog, "z", dd->nc[ZZ], options.cellSizeZ);
2548     }
2549
2550     /* Set the multi-body cut-off and cellsize limit for DLB */
2551     comm->cutoff_mbody   = systemInfo.minCutoffForMultiBody;
2552     comm->cellsize_limit = systemInfo.cellsizeLimit;
2553     if (systemInfo.haveInterDomainBondeds && systemInfo.increaseMultiBodyCutoff)
2554     {
2555         if (systemInfo.filterBondedCommunication || !isDlbDisabled(comm))
2556         {
2557             /* Set the bonded communication distance to halfway
2558              * the minimum and the maximum,
2559              * since the extra communication cost is nearly zero.
2560              */
2561             real acs           = average_cellsize_min(ddbox, dd->nc);
2562             comm->cutoff_mbody = 0.5*(systemInfo.minCutoffForMultiBody + acs);
2563             if (!isDlbDisabled(comm))
2564             {
2565                 /* Check if this does not limit the scaling */
2566                 comm->cutoff_mbody = std::min(comm->cutoff_mbody,
2567                                               options.dlbScaling*acs);
2568             }
2569             if (!systemInfo.filterBondedCommunication)
2570             {
2571                 /* Without bBondComm do not go beyond the n.b. cut-off */
2572                 comm->cutoff_mbody = std::min(comm->cutoff_mbody, systemInfo.cutoff);
2573                 if (comm->cellsize_limit >= systemInfo.cutoff)
2574                 {
2575                     /* We don't loose a lot of efficieny
2576                      * when increasing it to the n.b. cut-off.
2577                      * It can even be slightly faster, because we need
2578                      * less checks for the communication setup.
2579                      */
2580                     comm->cutoff_mbody = systemInfo.cutoff;
2581                 }
2582             }
2583             /* Check if we did not end up below our original limit */
2584             comm->cutoff_mbody = std::max(comm->cutoff_mbody,
2585                                           systemInfo.minCutoffForMultiBody);
2586
2587             if (comm->cutoff_mbody > comm->cellsize_limit)
2588             {
2589                 comm->cellsize_limit = comm->cutoff_mbody;
2590             }
2591         }
2592         /* Without DLB and cutoff_mbody<cutoff, cutoff_mbody is dynamic */
2593     }
2594
2595     if (debug)
2596     {
2597         fprintf(debug, "Bonded atom communication beyond the cut-off: %s\n"
2598                 "cellsize limit %f\n",
2599                 gmx::boolToString(systemInfo.filterBondedCommunication),
2600                 comm->cellsize_limit);
2601     }
2602
2603     if (MASTER(cr))
2604     {
2605         check_dd_restrictions(dd, ir, mdlog);
2606     }
2607 }
2608
2609 static char *init_bLocalCG(const gmx_mtop_t *mtop)
2610 {
2611     int   ncg, cg;
2612     char *bLocalCG;
2613
2614     ncg = ncg_mtop(mtop);
2615     snew(bLocalCG, ncg);
2616     for (cg = 0; cg < ncg; cg++)
2617     {
2618         bLocalCG[cg] = FALSE;
2619     }
2620
2621     return bLocalCG;
2622 }
2623
2624 void dd_init_bondeds(FILE *fplog,
2625                      gmx_domdec_t *dd,
2626                      const gmx_mtop_t *mtop,
2627                      const gmx_vsite_t *vsite,
2628                      const t_inputrec *ir,
2629                      gmx_bool bBCheck, cginfo_mb_t *cginfo_mb)
2630 {
2631     gmx_domdec_comm_t *comm;
2632
2633     dd_make_reverse_top(fplog, dd, mtop, vsite, ir, bBCheck);
2634
2635     comm = dd->comm;
2636
2637     if (comm->systemInfo.filterBondedCommunication)
2638     {
2639         /* Communicate atoms beyond the cut-off for bonded interactions */
2640         comm = dd->comm;
2641
2642         comm->cglink = make_charge_group_links(mtop, dd, cginfo_mb);
2643
2644         comm->bLocalCG = init_bLocalCG(mtop);
2645     }
2646     else
2647     {
2648         /* Only communicate atoms based on cut-off */
2649         comm->cglink   = nullptr;
2650         comm->bLocalCG = nullptr;
2651     }
2652 }
2653
2654 static void writeSettings(gmx::TextWriter       *log,
2655                           gmx_domdec_t          *dd,
2656                           const gmx_mtop_t      *mtop,
2657                           const t_inputrec      *ir,
2658                           gmx_bool               bDynLoadBal,
2659                           real                   dlb_scale,
2660                           const gmx_ddbox_t     *ddbox)
2661 {
2662     gmx_domdec_comm_t *comm;
2663     int                d;
2664     ivec               np;
2665     real               limit, shrink;
2666
2667     comm = dd->comm;
2668
2669     if (bDynLoadBal)
2670     {
2671         log->writeString("The maximum number of communication pulses is:");
2672         for (d = 0; d < dd->ndim; d++)
2673         {
2674             log->writeStringFormatted(" %c %d", dim2char(dd->dim[d]), comm->cd[d].np_dlb);
2675         }
2676         log->ensureLineBreak();
2677         log->writeLineFormatted("The minimum size for domain decomposition cells is %.3f nm", comm->cellsize_limit);
2678         log->writeLineFormatted("The requested allowed shrink of DD cells (option -dds) is: %.2f", dlb_scale);
2679         log->writeString("The allowed shrink of domain decomposition cells is:");
2680         for (d = 0; d < DIM; d++)
2681         {
2682             if (dd->nc[d] > 1)
2683             {
2684                 if (d >= ddbox->npbcdim && dd->nc[d] == 2)
2685                 {
2686                     shrink = 0;
2687                 }
2688                 else
2689                 {
2690                     shrink =
2691                         comm->cellsize_min_dlb[d]/
2692                         (ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
2693                 }
2694                 log->writeStringFormatted(" %c %.2f", dim2char(d), shrink);
2695             }
2696         }
2697         log->ensureLineBreak();
2698     }
2699     else
2700     {
2701         set_dd_cell_sizes_slb(dd, ddbox, setcellsizeslbPULSE_ONLY, np);
2702         log->writeString("The initial number of communication pulses is:");
2703         for (d = 0; d < dd->ndim; d++)
2704         {
2705             log->writeStringFormatted(" %c %d", dim2char(dd->dim[d]), np[dd->dim[d]]);
2706         }
2707         log->ensureLineBreak();
2708         log->writeString("The initial domain decomposition cell size is:");
2709         for (d = 0; d < DIM; d++)
2710         {
2711             if (dd->nc[d] > 1)
2712             {
2713                 log->writeStringFormatted(" %c %.2f nm",
2714                                           dim2char(d), dd->comm->cellsize_min[d]);
2715             }
2716         }
2717         log->ensureLineBreak();
2718         log->writeLine();
2719     }
2720
2721     const bool haveInterDomainVsites =
2722         (countInterUpdategroupVsites(*mtop, comm->systemInfo.updateGroupingPerMoleculetype) != 0);
2723
2724     if (comm->systemInfo.haveInterDomainBondeds ||
2725         haveInterDomainVsites ||
2726         comm->systemInfo.haveSplitConstraints ||
2727         comm->systemInfo.haveSplitSettles)
2728     {
2729         std::string decompUnits;
2730         if (comm->systemInfo.useUpdateGroups)
2731         {
2732             decompUnits = "atom groups";
2733         }
2734         else
2735         {
2736             decompUnits = "atoms";
2737         }
2738
2739         log->writeLineFormatted("The maximum allowed distance for %s involved in interactions is:", decompUnits.c_str());
2740         log->writeLineFormatted("%40s  %-7s %6.3f nm", "non-bonded interactions", "", comm->systemInfo.cutoff);
2741
2742         if (bDynLoadBal)
2743         {
2744             limit = dd->comm->cellsize_limit;
2745         }
2746         else
2747         {
2748             if (dd->unitCellInfo.ddBoxIsDynamic)
2749             {
2750                 log->writeLine("(the following are initial values, they could change due to box deformation)");
2751             }
2752             limit = dd->comm->cellsize_min[XX];
2753             for (d = 1; d < DIM; d++)
2754             {
2755                 limit = std::min(limit, dd->comm->cellsize_min[d]);
2756             }
2757         }
2758
2759         if (comm->systemInfo.haveInterDomainBondeds)
2760         {
2761             log->writeLineFormatted("%40s  %-7s %6.3f nm",
2762                                     "two-body bonded interactions", "(-rdd)",
2763                                     std::max(comm->systemInfo.cutoff, comm->cutoff_mbody));
2764             log->writeLineFormatted("%40s  %-7s %6.3f nm",
2765                                     "multi-body bonded interactions", "(-rdd)",
2766                                     (comm->systemInfo.filterBondedCommunication || isDlbOn(dd->comm)) ? comm->cutoff_mbody : std::min(comm->systemInfo.cutoff, limit));
2767         }
2768         if (haveInterDomainVsites)
2769         {
2770             log->writeLineFormatted("%40s  %-7s %6.3f nm",
2771                                     "virtual site constructions", "(-rcon)", limit);
2772         }
2773         if (comm->systemInfo.haveSplitConstraints || comm->systemInfo.haveSplitSettles)
2774         {
2775             std::string separation = gmx::formatString("atoms separated by up to %d constraints",
2776                                                        1+ir->nProjOrder);
2777             log->writeLineFormatted("%40s  %-7s %6.3f nm\n",
2778                                     separation.c_str(), "(-rcon)", limit);
2779         }
2780         log->ensureLineBreak();
2781     }
2782 }
2783
2784 static void logSettings(const gmx::MDLogger &mdlog,
2785                         gmx_domdec_t        *dd,
2786                         const gmx_mtop_t    *mtop,
2787                         const t_inputrec    *ir,
2788                         real                 dlb_scale,
2789                         const gmx_ddbox_t   *ddbox)
2790 {
2791     gmx::StringOutputStream stream;
2792     gmx::TextWriter         log(&stream);
2793     writeSettings(&log, dd, mtop, ir, isDlbOn(dd->comm), dlb_scale, ddbox);
2794     if (dd->comm->dlbState == DlbState::offCanTurnOn)
2795     {
2796         {
2797             log.ensureEmptyLine();
2798             log.writeLine("When dynamic load balancing gets turned on, these settings will change to:");
2799         }
2800         writeSettings(&log, dd, mtop, ir, true, dlb_scale, ddbox);
2801     }
2802     GMX_LOG(mdlog.info).asParagraph().appendText(stream.toString());
2803 }
2804
2805 static void set_cell_limits_dlb(const gmx::MDLogger &mdlog,
2806                                 gmx_domdec_t        *dd,
2807                                 real                 dlb_scale,
2808                                 const t_inputrec    *ir,
2809                                 const gmx_ddbox_t   *ddbox)
2810 {
2811     gmx_domdec_comm_t *comm;
2812     int                d, dim, npulse, npulse_d_max, npulse_d;
2813     gmx_bool           bNoCutOff;
2814
2815     comm = dd->comm;
2816
2817     bNoCutOff = (ir->rvdw == 0 || ir->rcoulomb == 0);
2818
2819     /* Determine the maximum number of comm. pulses in one dimension */
2820
2821     comm->cellsize_limit = std::max(comm->cellsize_limit, comm->cutoff_mbody);
2822
2823     /* Determine the maximum required number of grid pulses */
2824     if (comm->cellsize_limit >= comm->systemInfo.cutoff)
2825     {
2826         /* Only a single pulse is required */
2827         npulse = 1;
2828     }
2829     else if (!bNoCutOff && comm->cellsize_limit > 0)
2830     {
2831         /* We round down slightly here to avoid overhead due to the latency
2832          * of extra communication calls when the cut-off
2833          * would be only slightly longer than the cell size.
2834          * Later cellsize_limit is redetermined,
2835          * so we can not miss interactions due to this rounding.
2836          */
2837         npulse = static_cast<int>(0.96 + comm->systemInfo.cutoff/comm->cellsize_limit);
2838     }
2839     else
2840     {
2841         /* There is no cell size limit */
2842         npulse = std::max(dd->nc[XX]-1, std::max(dd->nc[YY]-1, dd->nc[ZZ]-1));
2843     }
2844
2845     if (!bNoCutOff && npulse > 1)
2846     {
2847         /* See if we can do with less pulses, based on dlb_scale */
2848         npulse_d_max = 0;
2849         for (d = 0; d < dd->ndim; d++)
2850         {
2851             dim      = dd->dim[d];
2852             npulse_d = static_cast<int>(1 + dd->nc[dim]*comm->systemInfo.cutoff
2853                                         /(ddbox->box_size[dim]*ddbox->skew_fac[dim]*dlb_scale));
2854             npulse_d_max = std::max(npulse_d_max, npulse_d);
2855         }
2856         npulse = std::min(npulse, npulse_d_max);
2857     }
2858
2859     /* This env var can override npulse */
2860     d = dd_getenv(mdlog, "GMX_DD_NPULSE", 0);
2861     if (d > 0)
2862     {
2863         npulse = d;
2864     }
2865
2866     comm->maxpulse       = 1;
2867     comm->bVacDLBNoLimit = (ir->ePBC == epbcNONE);
2868     for (d = 0; d < dd->ndim; d++)
2869     {
2870         comm->cd[d].np_dlb    = std::min(npulse, dd->nc[dd->dim[d]]-1);
2871         comm->maxpulse        = std::max(comm->maxpulse, comm->cd[d].np_dlb);
2872         if (comm->cd[d].np_dlb < dd->nc[dd->dim[d]]-1)
2873         {
2874             comm->bVacDLBNoLimit = FALSE;
2875         }
2876     }
2877
2878     /* cellsize_limit is set for LINCS in init_domain_decomposition */
2879     if (!comm->bVacDLBNoLimit)
2880     {
2881         comm->cellsize_limit = std::max(comm->cellsize_limit,
2882                                         comm->systemInfo.cutoff/comm->maxpulse);
2883     }
2884     comm->cellsize_limit = std::max(comm->cellsize_limit, comm->cutoff_mbody);
2885     /* Set the minimum cell size for each DD dimension */
2886     for (d = 0; d < dd->ndim; d++)
2887     {
2888         if (comm->bVacDLBNoLimit ||
2889             comm->cd[d].np_dlb*comm->cellsize_limit >= comm->systemInfo.cutoff)
2890         {
2891             comm->cellsize_min_dlb[dd->dim[d]] = comm->cellsize_limit;
2892         }
2893         else
2894         {
2895             comm->cellsize_min_dlb[dd->dim[d]] =
2896                 comm->systemInfo.cutoff/comm->cd[d].np_dlb;
2897         }
2898     }
2899     if (comm->cutoff_mbody <= 0)
2900     {
2901         comm->cutoff_mbody = std::min(comm->systemInfo.cutoff, comm->cellsize_limit);
2902     }
2903     if (isDlbOn(comm))
2904     {
2905         set_dlb_limits(dd);
2906     }
2907 }
2908
2909 gmx_bool dd_bonded_molpbc(const gmx_domdec_t *dd, int ePBC)
2910 {
2911     /* If each molecule is a single charge group
2912      * or we use domain decomposition for each periodic dimension,
2913      * we do not need to take pbc into account for the bonded interactions.
2914      */
2915     return (ePBC != epbcNONE && dd->comm->systemInfo.haveInterDomainBondeds &&
2916             !(dd->nc[XX] > 1 &&
2917               dd->nc[YY] > 1 &&
2918               (dd->nc[ZZ] > 1 || ePBC == epbcXY)));
2919 }
2920
2921 /*! \brief Sets grid size limits and PP-PME setup, prints settings to log */
2922 static void set_ddgrid_parameters(const gmx::MDLogger &mdlog,
2923                                   gmx_domdec_t *dd, real dlb_scale,
2924                                   const gmx_mtop_t *mtop, const t_inputrec *ir,
2925                                   const gmx_ddbox_t *ddbox)
2926 {
2927     gmx_domdec_comm_t *comm        = dd->comm;
2928     DDRankSetup       &ddRankSetup = comm->ddRankSetup;
2929
2930     if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
2931     {
2932         init_ddpme(dd, &ddRankSetup.ddpme[0], 0);
2933         if (ddRankSetup.npmedecompdim >= 2)
2934         {
2935             init_ddpme(dd, &ddRankSetup.ddpme[1], 1);
2936         }
2937     }
2938     else
2939     {
2940         ddRankSetup.npmenodes = 0;
2941         if (dd->pme_nodeid >= 0)
2942         {
2943             gmx_fatal_collective(FARGS, dd->mpi_comm_all, DDMASTER(dd),
2944                                  "Can not have separate PME ranks without PME electrostatics");
2945         }
2946     }
2947
2948     if (debug)
2949     {
2950         fprintf(debug, "The DD cut-off is %f\n", comm->systemInfo.cutoff);
2951     }
2952     if (!isDlbDisabled(comm))
2953     {
2954         set_cell_limits_dlb(mdlog, dd, dlb_scale, ir, ddbox);
2955     }
2956
2957     logSettings(mdlog, dd, mtop, ir, dlb_scale, ddbox);
2958
2959     real vol_frac;
2960     if (ir->ePBC == epbcNONE)
2961     {
2962         vol_frac = 1 - 1/static_cast<double>(dd->nnodes);
2963     }
2964     else
2965     {
2966         vol_frac =
2967             (1 + comm_box_frac(dd->nc, comm->systemInfo.cutoff, ddbox))/static_cast<double>(dd->nnodes);
2968     }
2969     if (debug)
2970     {
2971         fprintf(debug, "Volume fraction for all DD zones: %f\n", vol_frac);
2972     }
2973     int natoms_tot = comm->cgs_gl.index[comm->cgs_gl.nr];
2974
2975     dd->ga2la      = new gmx_ga2la_t(natoms_tot,
2976                                      static_cast<int>(vol_frac*natoms_tot));
2977 }
2978
2979 /*! \brief Get some important DD parameters which can be modified by env.vars */
2980 static DDSettings
2981 getDDSettings(const gmx::MDLogger     &mdlog,
2982               const DomdecOptions     &options,
2983               const gmx::MdrunOptions &mdrunOptions,
2984               const t_inputrec        &ir)
2985 {
2986     DDSettings ddSettings;
2987
2988     ddSettings.useSendRecv2        = (dd_getenv(mdlog, "GMX_DD_USE_SENDRECV2", 0) != 0);
2989     ddSettings.dlb_scale_lim       = dd_getenv(mdlog, "GMX_DLB_MAX_BOX_SCALING", 10);
2990     ddSettings.useDDOrderZYX       = bool(dd_getenv(mdlog, "GMX_DD_ORDER_ZYX", 0));
2991     ddSettings.useCartesianReorder = bool(dd_getenv(mdlog, "GMX_NO_CART_REORDER", 1));
2992     ddSettings.eFlop               = dd_getenv(mdlog, "GMX_DLB_BASED_ON_FLOPS", 0);
2993     const int recload              = dd_getenv(mdlog, "GMX_DD_RECORD_LOAD", 1);
2994     ddSettings.nstDDDump           = dd_getenv(mdlog, "GMX_DD_NST_DUMP", 0);
2995     ddSettings.nstDDDumpGrid       = dd_getenv(mdlog, "GMX_DD_NST_DUMP_GRID", 0);
2996     ddSettings.DD_debug            = dd_getenv(mdlog, "GMX_DD_DEBUG", 0);
2997
2998     if (ddSettings.useSendRecv2)
2999     {
3000         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");
3001     }
3002
3003     if (ddSettings.eFlop)
3004     {
3005         GMX_LOG(mdlog.info).appendText("Will load balance based on FLOP count");
3006         ddSettings.recordLoad = true;
3007     }
3008     else
3009     {
3010         ddSettings.recordLoad = (wallcycle_have_counter() && recload > 0);
3011     }
3012
3013     ddSettings.initialDlbState =
3014         determineInitialDlbState(mdlog, options.dlbOption, ddSettings.recordLoad, mdrunOptions, &ir);
3015     GMX_LOG(mdlog.info).appendTextFormatted("Dynamic load balancing: %s",
3016                                             edlbs_names[static_cast<int>(ddSettings.initialDlbState)]);
3017
3018     return ddSettings;
3019 }
3020
3021 gmx_domdec_t::gmx_domdec_t(const t_inputrec &ir) :
3022     unitCellInfo(ir)
3023 {
3024 }
3025
3026 gmx_domdec_t *init_domain_decomposition(const gmx::MDLogger           &mdlog,
3027                                         t_commrec                     *cr,
3028                                         const DomdecOptions           &options,
3029                                         const gmx::MdrunOptions       &mdrunOptions,
3030                                         const gmx_mtop_t              *mtop,
3031                                         const t_inputrec              *ir,
3032                                         const matrix                   box,
3033                                         gmx::ArrayRef<const gmx::RVec> xGlobal,
3034                                         gmx::LocalAtomSetManager      *atomSets)
3035 {
3036     GMX_LOG(mdlog.info).appendTextFormatted(
3037             "\nInitializing Domain Decomposition on %d ranks", cr->nnodes);
3038
3039     DDSettings  ddSettings = getDDSettings(mdlog, options, mdrunOptions, *ir);
3040     if (ddSettings.eFlop > 1)
3041     {
3042         /* Ensure that we have different random flop counts on different ranks */
3043         srand(1 + cr->nodeid);
3044     }
3045
3046     DDSystemInfo systemInfo = getSystemInfo(mdlog, cr, options, mtop, ir, box, xGlobal);
3047
3048     gmx_ddbox_t  ddbox      = {0};
3049     DDSetup      ddSetup    = getDDSetup(mdlog, cr, options, ddSettings, systemInfo,
3050                                          mtop, ir, box, xGlobal, &ddbox);
3051
3052     cr->npmenodes = ddSetup.numPmeRanks;
3053
3054     DDRankSetup ddRankSetup = getDDRankSetup(mdlog, cr, ddSetup, *ir);
3055
3056     ivec        ddCellIndex = { 0, 0, 0 };
3057     makeGroupCommunicators(mdlog, ddSettings, ddSetup, options.rankOrder,
3058                            &ddRankSetup, cr, ddCellIndex);
3059
3060     gmx_domdec_t *dd = new gmx_domdec_t(*ir);
3061
3062     copy_ivec(ddCellIndex, dd->ci);
3063
3064     dd->comm = init_dd_comm();
3065
3066     dd->comm->ddRankSetup = ddRankSetup;
3067
3068     set_dd_limits(mdlog, cr, dd, options,
3069                   ddSettings, systemInfo, ddSetup,
3070                   mtop, ir,
3071                   ddbox);
3072
3073     setupGroupCommunication(mdlog, ddSettings, cr, dd);
3074
3075     if (thisRankHasDuty(cr, DUTY_PP))
3076     {
3077         set_ddgrid_parameters(mdlog, dd, options.dlbScaling, mtop, ir, &ddbox);
3078
3079         setup_neighbor_relations(dd);
3080     }
3081
3082     /* Set overallocation to avoid frequent reallocation of arrays */
3083     set_over_alloc_dd(TRUE);
3084
3085     dd->atomSets = atomSets;
3086
3087     return dd;
3088 }
3089
3090 static gmx_bool test_dd_cutoff(t_commrec                     *cr,
3091                                const matrix                   box,
3092                                gmx::ArrayRef<const gmx::RVec> x,
3093                                real                           cutoffRequested)
3094 {
3095     gmx_domdec_t *dd;
3096     gmx_ddbox_t   ddbox;
3097     int           d, dim, np;
3098     real          inv_cell_size;
3099     int           LocallyLimited;
3100
3101     dd = cr->dd;
3102
3103     set_ddbox(*dd, false, box, true, x, &ddbox);
3104
3105     LocallyLimited = 0;
3106
3107     for (d = 0; d < dd->ndim; d++)
3108     {
3109         dim = dd->dim[d];
3110
3111         inv_cell_size = DD_CELL_MARGIN*dd->nc[dim]/ddbox.box_size[dim];
3112         if (dd->unitCellInfo.ddBoxIsDynamic)
3113         {
3114             inv_cell_size *= DD_PRES_SCALE_MARGIN;
3115         }
3116
3117         np = 1 + static_cast<int>(cutoffRequested*inv_cell_size*ddbox.skew_fac[dim]);
3118
3119         if (!isDlbDisabled(dd->comm) && (dim < ddbox.npbcdim) && (dd->comm->cd[d].np_dlb > 0))
3120         {
3121             if (np > dd->comm->cd[d].np_dlb)
3122             {
3123                 return FALSE;
3124             }
3125
3126             /* If a current local cell size is smaller than the requested
3127              * cut-off, we could still fix it, but this gets very complicated.
3128              * Without fixing here, we might actually need more checks.
3129              */
3130             real cellSizeAlongDim = (dd->comm->cell_x1[dim] - dd->comm->cell_x0[dim])*ddbox.skew_fac[dim];
3131             if (cellSizeAlongDim*dd->comm->cd[d].np_dlb < cutoffRequested)
3132             {
3133                 LocallyLimited = 1;
3134             }
3135         }
3136     }
3137
3138     if (!isDlbDisabled(dd->comm))
3139     {
3140         /* If DLB is not active yet, we don't need to check the grid jumps.
3141          * Actually we shouldn't, because then the grid jump data is not set.
3142          */
3143         if (isDlbOn(dd->comm) &&
3144             check_grid_jump(0, dd, cutoffRequested, &ddbox, FALSE))
3145         {
3146             LocallyLimited = 1;
3147         }
3148
3149         gmx_sumi(1, &LocallyLimited, cr);
3150
3151         if (LocallyLimited > 0)
3152         {
3153             return FALSE;
3154         }
3155     }
3156
3157     return TRUE;
3158 }
3159
3160 gmx_bool change_dd_cutoff(t_commrec                     *cr,
3161                           const matrix                   box,
3162                           gmx::ArrayRef<const gmx::RVec> x,
3163                           real                           cutoffRequested)
3164 {
3165     gmx_bool bCutoffAllowed;
3166
3167     bCutoffAllowed = test_dd_cutoff(cr, box, x, cutoffRequested);
3168
3169     if (bCutoffAllowed)
3170     {
3171         cr->dd->comm->systemInfo.cutoff = cutoffRequested;
3172     }
3173
3174     return bCutoffAllowed;
3175 }