src/gromacs/domdec/domdec.cpp

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