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