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