src/gromacs/domdec/domdec_setup.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2008,2009,2010,2011,2012 by the GROMACS development team.
   5  * Copyright (c) 2013,2014,2015,2017,2018 by the GROMACS development team.
   6  * Copyright (c) 2019,2020,2021, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37
  38 /*! \internal \file
  39  *
  40  * \brief This file defines functions used by the domdec module
  41  * in its initial setup phase.
  42  *
  43  * \author Berk Hess <hess@kth.se>
  44  * \ingroup module_domdec
  45  */
  46
  47 #include "gmxpre.h"
  48
  49 #include "domdec_setup.h"
  50
  51 #include <cassert>
  52 #include <cinttypes>
  53 #include <cmath>
  54 #include <cstdio>
  55
  56 #include <numeric>
  57
  58 #include "gromacs/domdec/domdec.h"
  59 #include "gromacs/domdec/domdec_struct.h"
  60 #include "gromacs/domdec/options.h"
  61 #include "gromacs/ewald/pme.h"
  62 #include "gromacs/gmxlib/network.h"
  63 #include "gromacs/math/units.h"
  64 #include "gromacs/math/utilities.h"
  65 #include "gromacs/math/vec.h"
  66 #include "gromacs/mdlib/perf_est.h"
  67 #include "gromacs/mdtypes/commrec.h"
  68 #include "gromacs/mdtypes/inputrec.h"
  69 #include "gromacs/mdtypes/md_enums.h"
  70 #include "gromacs/pbcutil/pbc.h"
  71 #include "gromacs/topology/topology.h"
  72 #include "gromacs/utility/fatalerror.h"
  73 #include "gromacs/utility/logger.h"
  74 #include "gromacs/utility/stringutil.h"
  75
  76 #include "box.h"
  77 #include "domdec_internal.h"
  78 #include "utility.h"
  79
  80 // TODO remove this when moving domdec into gmx namespace
  81 using gmx::DomdecOptions;
  82
  83 /*! \brief Margin for setting up the DD grid */
  84 #define DD_GRID_MARGIN_PRES_SCALE 1.05
  85
  86 /*! \brief Factorize \p n.
  87  *
  88  * \param[in]    n     Value to factorize
  89  * \param[out]   fac   Vector of factors (to be allocated in this function)
  90  * \param[out]   mfac  Vector with the number of times each factor repeats in the factorization (to be allocated in this function)
  91  */
  92 static void factorize(int n, std::vector<int>* fac, std::vector<int>* mfac)
  93 {
  94     if (n <= 0)
  95     {
  96         gmx_fatal(FARGS, "Can only factorize positive integers.");
  97     }
  98
  99     /* Decompose n in factors */
 100     fac->clear();
 101     mfac->clear();
 102     int d = 2;
 103     while (n > 1)
 104     {
 105         while (n % d == 0)
 106         {
 107             if (fac->empty() || fac->back() != d)
 108             {
 109                 fac->push_back(d);
 110                 mfac->push_back(1);
 111             }
 112             else
 113             {
 114                 mfac->back()++;
 115             }
 116             n /= d;
 117         }
 118         d++;
 119     }
 120 }
 121
 122 /*! \brief Find largest divisor of \p n smaller than \p n*/
 123 static int largest_divisor(int n)
 124 {
 125     std::vector<int> div;
 126     std::vector<int> mdiv;
 127     factorize(n, &div, &mdiv);
 128
 129     return div.back();
 130 }
 131
 132 /*! \brief Returns TRUE when there are enough PME ranks for the ratio */
 133 static gmx_bool fits_pme_ratio(int nrank_tot, int nrank_pme, float ratio)
 134 {
 135     return (static_cast<double>(nrank_pme) / static_cast<double>(nrank_tot) > 0.95 * ratio);
 136 }
 137
 138 /*! \brief Returns TRUE when npme out of ntot ranks doing PME is expected to give reasonable performance */
 139 static gmx_bool fits_pp_pme_perf(int ntot, int npme, float ratio)
 140 {
 141     std::vector<int> div;
 142     std::vector<int> mdiv;
 143     factorize(ntot - npme, &div, &mdiv);
 144
 145     int npp_root3  = gmx::roundToInt(std::cbrt(ntot - npme));
 146     int npme_root2 = gmx::roundToInt(std::sqrt(static_cast<double>(npme)));
 147
 148     /* The check below gives a reasonable division:
 149      * factor 5 allowed at 5 or more PP ranks,
 150      * factor 7 allowed at 49 or more PP ranks.
 151      */
 152     if (div.back() > 3 + npp_root3)
 153     {
 154         return FALSE;
 155     }
 156
 157     /* Check if the number of PP and PME ranks have a reasonable sized
 158      * denominator in common, such that we can use 2D PME decomposition
 159      * when required (which requires nx_pp == nx_pme).
 160      * The factor of 2 allows for a maximum ratio of 2^2=4
 161      * between nx_pme and ny_pme.
 162      */
 163     if (std::gcd(ntot - npme, npme) * 2 < npme_root2)
 164     {
 165         return FALSE;
 166     }
 167
 168     /* Does this division gives a reasonable PME load? */
 169     return fits_pme_ratio(ntot, npme, ratio);
 170 }
 171
 172 /*! \brief Make a guess for the number of PME ranks to use. */
 173 static int guess_npme(const gmx::MDLogger& mdlog,
 174                       const gmx_mtop_t&    mtop,
 175                       const t_inputrec&    ir,
 176                       const matrix         box,
 177                       int                  nrank_tot)
 178 {
 179     float ratio = pme_load_estimate(mtop, ir, box);
 180
 181     GMX_LOG(mdlog.info).appendTextFormatted("Guess for relative PME load: %.2f", ratio);
 182
 183     /* We assume the optimal rank ratio is close to the load ratio.
 184      * The communication load is neglected,
 185      * but (hopefully) this will balance out between PP and PME.
 186      */
 187
 188     if (!fits_pme_ratio(nrank_tot, nrank_tot / 2, ratio))
 189     {
 190         /* We would need more than nrank_tot/2 PME only nodes,
 191          * which is not possible. Since the PME load is very high,
 192          * we will not loose much performance when all ranks do PME.
 193          */
 194
 195         return 0;
 196     }
 197
 198     /* First try to find npme as a factor of nrank_tot up to nrank_tot/3.
 199      * We start with a minimum PME node fraction of 1/16
 200      * and avoid ratios which lead to large prime factors in nnodes-npme.
 201      */
 202     int npme = (nrank_tot + 15) / 16;
 203     while (npme <= nrank_tot / 3)
 204     {
 205         if (nrank_tot % npme == 0)
 206         {
 207             /* Note that fits_perf might change the PME grid,
 208              * in the current implementation it does not.
 209              */
 210             if (fits_pp_pme_perf(nrank_tot, npme, ratio))
 211             {
 212                 break;
 213             }
 214         }
 215         npme++;
 216     }
 217     if (npme > nrank_tot / 3)
 218     {
 219         /* Try any possible number for npme */
 220         npme = 1;
 221         while (npme <= nrank_tot / 2)
 222         {
 223             /* Note that fits_perf may change the PME grid */
 224             if (fits_pp_pme_perf(nrank_tot, npme, ratio))
 225             {
 226                 break;
 227             }
 228             npme++;
 229         }
 230     }
 231     if (npme > nrank_tot / 2)
 232     {
 233         gmx_fatal(FARGS,
 234                   "Could not find an appropriate number of separate PME ranks. i.e. >= %5f*#ranks "
 235                   "(%d) and <= #ranks/2 (%d) and reasonable performance wise (grid_x=%d, "
 236                   "grid_y=%d).\n"
 237                   "Use the -npme option of mdrun or change the number of ranks or the PME grid "
 238                   "dimensions, see the manual for details.",
 239                   ratio,
 240                   gmx::roundToInt(0.95 * ratio * nrank_tot),
 241                   nrank_tot / 2,
 242                   ir.nkx,
 243                   ir.nky);
 244     }
 245     else
 246     {
 247         GMX_LOG(mdlog.info)
 248                 .appendTextFormatted(
 249                         "Will use %d particle-particle and %d PME only ranks\n"
 250                         "This is a guess, check the performance at the end of the log file",
 251                         nrank_tot - npme,
 252                         npme);
 253     }
 254
 255     return npme;
 256 }
 257
 258 /*! \brief Return \p n divided by \p f rounded up to the next integer. */
 259 static int div_up(int n, int f)
 260 {
 261     return (n + f - 1) / f;
 262 }
 263
 264 real comm_box_frac(const gmx::IVec& dd_nc, real cutoff, const gmx_ddbox_t& ddbox)
 265 {
 266     rvec nw;
 267
 268     for (int i = 0; i < DIM; i++)
 269     {
 270         real bt = ddbox.box_size[i] * ddbox.skew_fac[i];
 271         nw[i]   = dd_nc[i] * cutoff / bt;
 272     }
 273
 274     real comm_vol = 0;
 275     for (int i = 0; i < DIM; i++)
 276     {
 277         if (dd_nc[i] > 1)
 278         {
 279             comm_vol += nw[i];
 280             for (int j = i + 1; j < DIM; j++)
 281             {
 282                 if (dd_nc[j] > 1)
 283                 {
 284                     comm_vol += nw[i] * nw[j] * M_PI / 4;
 285                     for (int k = j + 1; k < DIM; k++)
 286                     {
 287                         if (dd_nc[k] > 1)
 288                         {
 289                             comm_vol += nw[i] * nw[j] * nw[k] * M_PI / 6;
 290                         }
 291                     }
 292                 }
 293             }
 294         }
 295     }
 296
 297     return comm_vol;
 298 }
 299
 300 /*! \brief Return whether the DD inhomogeneous in the z direction */
 301 static gmx_bool inhomogeneous_z(const t_inputrec& ir)
 302 {
 303     return ((EEL_PME(ir.coulombtype) || ir.coulombtype == CoulombInteractionType::Ewald)
 304             && ir.pbcType == PbcType::Xyz && ir.ewald_geometry == EwaldGeometry::ThreeDC);
 305 }
 306
 307 /*! \brief Estimate cost of PME FFT communication
 308  *
 309  * This only takes the communication into account and not imbalance
 310  * in the calculation. But the imbalance in communication and calculation
 311  * are similar and therefore these formulas also prefer load balance
 312  * in the FFT and pme_solve calculation.
 313  */
 314 static float comm_pme_cost_vol(int npme, int a, int b, int c)
 315 {
 316     /* We use a float here, since an integer might overflow */
 317     float comm_vol = npme - 1;
 318     comm_vol *= npme;
 319     comm_vol *= div_up(a, npme);
 320     comm_vol *= div_up(b, npme);
 321     comm_vol *= c;
 322
 323     return comm_vol;
 324 }
 325
 326 /*! \brief Estimate cost of communication for a possible domain decomposition. */
 327 static float comm_cost_est(real               limit,
 328                            real               cutoff,
 329                            const matrix       box,
 330                            const gmx_ddbox_t& ddbox,
 331                            int                natoms,
 332                            const t_inputrec&  ir,
 333                            float              pbcdxr,
 334                            int                npme_tot,
 335                            const gmx::IVec&   nc)
 336 {
 337     gmx::IVec npme = { 1, 1, 1 };
 338     rvec      bt;
 339     /* This is the cost of a pbc_dx call relative to the cost
 340      * of communicating the coordinate and force of an atom.
 341      * This will be machine dependent.
 342      * These factors are for x86 with SMP or Infiniband.
 343      */
 344     float pbcdx_rect_fac = 0.1;
 345     float pbcdx_tric_fac = 0.2;
 346
 347     /* Check the DD algorithm restrictions */
 348     if ((ir.pbcType == PbcType::XY && ir.nwall < 2 && nc[ZZ] > 1)
 349         || (ir.pbcType == PbcType::Screw && (nc[XX] == 1 || nc[YY] > 1 || nc[ZZ] > 1)))
 350     {
 351         return -1;
 352     }
 353
 354     if (inhomogeneous_z(ir) && nc[ZZ] > 1)
 355     {
 356         return -1;
 357     }
 358
 359     assert(ddbox.npbcdim <= DIM);
 360
 361     /* Check if the triclinic requirements are met */
 362     for (int i = 0; i < DIM; i++)
 363     {
 364         for (int j = i + 1; j < ddbox.npbcdim; j++)
 365         {
 366             if (box[j][i] != 0 || ir.deform[j][i] != 0
 367                 || (ir.epc != PressureCoupling::No && ir.compress[j][i] != 0))
 368             {
 369                 if (nc[j] > 1 && nc[i] == 1)
 370                 {
 371                     return -1;
 372                 }
 373             }
 374         }
 375     }
 376
 377     for (int i = 0; i < DIM; i++)
 378     {
 379         bt[i] = ddbox.box_size[i] * ddbox.skew_fac[i];
 380
 381         /* Without PBC and with 2 cells, there are no lower limits on the cell size */
 382         if (!(i >= ddbox.npbcdim && nc[i] <= 2) && bt[i] < nc[i] * limit)
 383         {
 384             return -1;
 385         }
 386         /* With PBC, check if the cut-off fits in nc[i]-1 cells */
 387         if (i < ddbox.npbcdim && nc[i] > 1 && (nc[i] - 1) * bt[i] < nc[i] * cutoff)
 388         {
 389             return -1;
 390         }
 391     }
 392
 393     if (npme_tot > 1)
 394     {
 395         /* The following choices should match those
 396          * in init_domain_decomposition in domdec.c.
 397          */
 398         if (nc[XX] == 1 && nc[YY] > 1)
 399         {
 400             npme[XX] = 1;
 401             npme[YY] = npme_tot;
 402         }
 403         else if (nc[YY] == 1)
 404         {
 405             npme[XX] = npme_tot;
 406             npme[YY] = 1;
 407         }
 408         else
 409         {
 410             /* Will we use 1D or 2D PME decomposition? */
 411             npme[XX] = (npme_tot % nc[XX] == 0) ? nc[XX] : npme_tot;
 412             npme[YY] = npme_tot / npme[XX];
 413         }
 414     }
 415
 416     if (EEL_PME(ir.coulombtype) || EVDW_PME(ir.vdwtype))
 417     {
 418         /* Check the PME grid restrictions.
 419          * Currently these can only be invalid here with too few grid lines
 420          * along the x dimension per rank doing PME.
 421          */
 422         int npme_x = (npme_tot > 1 ? npme[XX] : nc[XX]);
 423
 424         /* Currently we don't have the OpenMP thread count available here.
 425          * But with threads we have only tighter restrictions and it's
 426          * probably better anyhow to avoid settings where we need to reduce
 427          * grid lines over multiple ranks, as the thread check will do.
 428          */
 429         bool useThreads     = true;
 430         bool errorsAreFatal = false;
 431         if (!gmx_pme_check_restrictions(
 432                     ir.pme_order, ir.nkx, ir.nky, ir.nkz, npme_x, useThreads, errorsAreFatal))
 433         {
 434             return -1;
 435         }
 436     }
 437
 438     /* When two dimensions are (nearly) equal, use more cells
 439      * for the smallest index, so the decomposition does not
 440      * depend sensitively on the rounding of the box elements.
 441      */
 442     for (int i = 0; i < DIM; i++)
 443     {
 444         for (int j = i + 1; j < DIM; j++)
 445         {
 446             /* Check if the box size is nearly identical,
 447              * in that case we prefer nx > ny  and ny > nz.
 448              */
 449             if (std::fabs(bt[j] - bt[i]) < 0.01 * bt[i] && nc[j] > nc[i])
 450             {
 451                 /* The XX/YY check is a bit compact. If nc[YY]==npme[YY]
 452                  * this means the swapped nc has nc[XX]==npme[XX],
 453                  * and we can also swap X and Y for PME.
 454                  */
 455                 /* Check if dimension i and j are equivalent for PME.
 456                  * For x/y: if nc[YY]!=npme[YY], we can not swap x/y
 457                  * For y/z: we can not have PME decomposition in z
 458                  */
 459                 if (npme_tot <= 1
 460                     || !((i == XX && j == YY && nc[YY] != npme[YY]) || (i == YY && j == ZZ && npme[YY] > 1)))
 461                 {
 462                     return -1;
 463                 }
 464             }
 465         }
 466     }
 467
 468     /* This function determines only half of the communication cost.
 469      * All PP, PME and PP-PME communication is symmetric
 470      * and the "back"-communication cost is identical to the forward cost.
 471      */
 472
 473     float comm_vol = comm_box_frac(nc, cutoff, ddbox);
 474
 475     float comm_pme = 0;
 476     for (int i = 0; i < 2; i++)
 477     {
 478         /* Determine the largest volume for PME x/f redistribution */
 479         if (nc[i] % npme[i] != 0)
 480         {
 481             float comm_vol_xf =
 482                     (nc[i] > npme[i]) ? (npme[i] == 2 ? 1.0 / 3.0 : 0.5)
 483                                       : (1.0 - std::gcd(nc[i], npme[i]) / static_cast<double>(npme[i]));
 484             comm_pme += 3 * natoms * comm_vol_xf;
 485         }
 486
 487         /* Grid overlap communication */
 488         if (npme[i] > 1)
 489         {
 490             const int nk      = (i == 0 ? ir.nkx : ir.nky);
 491             const int overlap = (nk % npme[i] == 0 ? ir.pme_order - 1 : ir.pme_order);
 492             float     temp    = npme[i];
 493             temp *= overlap;
 494             temp *= ir.nkx;
 495             temp *= ir.nky;
 496             temp *= ir.nkz;
 497             temp /= nk;
 498             comm_pme += temp;
 499             /* Old line comm_pme += npme[i]*overlap*ir.nkx*ir.nky*ir.nkz/nk; */
 500         }
 501     }
 502
 503     comm_pme += comm_pme_cost_vol(npme[YY], ir.nky, ir.nkz, ir.nkx);
 504     comm_pme += comm_pme_cost_vol(npme[XX], ir.nkx, ir.nky, ir.nkz);
 505
 506     /* Add cost of pbc_dx for bondeds */
 507     float cost_pbcdx = 0;
 508     if ((nc[XX] == 1 || nc[YY] == 1) || (nc[ZZ] == 1 && ir.pbcType != PbcType::XY))
 509     {
 510         if ((ddbox.tric_dir[XX] && nc[XX] == 1) || (ddbox.tric_dir[YY] && nc[YY] == 1))
 511         {
 512             cost_pbcdx = pbcdxr * pbcdx_tric_fac;
 513         }
 514         else
 515         {
 516             cost_pbcdx = pbcdxr * pbcdx_rect_fac;
 517         }
 518     }
 519
 520     if (debug)
 521     {
 522         fprintf(debug,
 523                 "nc %2d %2d %2d %2d %2d vol pp %6.4f pbcdx %6.4f pme %9.3e tot %9.3e\n",
 524                 nc[XX],
 525                 nc[YY],
 526                 nc[ZZ],
 527                 npme[XX],
 528                 npme[YY],
 529                 comm_vol,
 530                 cost_pbcdx,
 531                 comm_pme / (3 * natoms),
 532                 comm_vol + cost_pbcdx + comm_pme / (3 * natoms));
 533     }
 534
 535     return 3 * natoms * (comm_vol + cost_pbcdx) + comm_pme;
 536 }
 537
 538 /*! \brief Assign penalty factors to possible domain decompositions,
 539  * based on the estimated communication costs. */
 540 static void assign_factors(const real         limit,
 541                            const real         cutoff,
 542                            const matrix       box,
 543                            const gmx_ddbox_t& ddbox,
 544                            int                natoms,
 545                            const t_inputrec&  ir,
 546                            float              pbcdxr,
 547                            int                npme,
 548                            int                ndiv,
 549                            const int*         div,
 550                            const int*         mdiv,
 551                            gmx::IVec*         irTryPtr,
 552                            gmx::IVec*         opt)
 553 {
 554     gmx::IVec& ir_try = *irTryPtr;
 555
 556     if (ndiv == 0)
 557     {
 558
 559         const float ce = comm_cost_est(limit, cutoff, box, ddbox, natoms, ir, pbcdxr, npme, ir_try);
 560         if (ce >= 0
 561             && ((*opt)[XX] == 0
 562                 || ce < comm_cost_est(limit, cutoff, box, ddbox, natoms, ir, pbcdxr, npme, *opt)))
 563         {
 564             *opt = ir_try;
 565         }
 566
 567         return;
 568     }
 569
 570     for (int x = mdiv[0]; x >= 0; x--)
 571     {
 572         for (int i = 0; i < x; i++)
 573         {
 574             ir_try[XX] *= div[0];
 575         }
 576         for (int y = mdiv[0] - x; y >= 0; y--)
 577         {
 578             for (int i = 0; i < y; i++)
 579             {
 580                 ir_try[YY] *= div[0];
 581             }
 582             for (int i = 0; i < mdiv[0] - x - y; i++)
 583             {
 584                 ir_try[ZZ] *= div[0];
 585             }
 586
 587             /* recurse */
 588             assign_factors(
 589                     limit, cutoff, box, ddbox, natoms, ir, pbcdxr, npme, ndiv - 1, div + 1, mdiv + 1, irTryPtr, opt);
 590
 591             for (int i = 0; i < mdiv[0] - x - y; i++)
 592             {
 593                 ir_try[ZZ] /= div[0];
 594             }
 595             for (int i = 0; i < y; i++)
 596             {
 597                 ir_try[YY] /= div[0];
 598             }
 599         }
 600         for (int i = 0; i < x; i++)
 601         {
 602             ir_try[XX] /= div[0];
 603         }
 604     }
 605 }
 606
 607 /*! \brief Determine the optimal distribution of DD cells for the
 608  * simulation system and number of MPI ranks
 609  *
 610  * \returns The optimal grid cell choice. The latter will contain all
 611  *          zeros if no valid cell choice exists. */
 612 static gmx::IVec optimizeDDCells(const gmx::MDLogger& mdlog,
 613                                  const int            numRanksRequested,
 614                                  const int            numPmeOnlyRanks,
 615                                  const real           cellSizeLimit,
 616                                  const gmx_mtop_t&    mtop,
 617                                  const matrix         box,
 618                                  const gmx_ddbox_t&   ddbox,
 619                                  const t_inputrec&    ir,
 620                                  const DDSystemInfo&  systemInfo)
 621 {
 622     double pbcdxr = 0;
 623
 624     const int numPPRanks = numRanksRequested - numPmeOnlyRanks;
 625
 626     GMX_LOG(mdlog.info)
 627             .appendTextFormatted(
 628                     "Optimizing the DD grid for %d cells with a minimum initial size of %.3f nm",
 629                     numPPRanks,
 630                     cellSizeLimit);
 631     if (inhomogeneous_z(ir))
 632     {
 633         GMX_LOG(mdlog.info)
 634                 .appendTextFormatted(
 635                         "Ewald_geometry=%s: assuming inhomogeneous particle distribution in z, "
 636                         "will not decompose in z.",
 637                         enumValueToString(ir.ewald_geometry));
 638     }
 639
 640
 641     // For cost estimates, we need the number of ranks doing PME work,
 642     // which is the number of PP ranks when not using separate
 643     // PME-only ranks.
 644     const int numRanksDoingPmeWork =
 645             (EEL_PME(ir.coulombtype) ? ((numPmeOnlyRanks > 0) ? numPmeOnlyRanks : numPPRanks) : 0);
 646
 647     if (systemInfo.haveInterDomainBondeds)
 648     {
 649         /* If we can skip PBC for distance calculations in plain-C bondeds,
 650          * we can save some time (e.g. 3D DD with pbc=xyz).
 651          * Here we ignore SIMD bondeds as they always do (fast) PBC.
 652          */
 653         count_bonded_distances(mtop, ir, &pbcdxr, nullptr);
 654         pbcdxr /= static_cast<double>(mtop.natoms);
 655     }
 656     else
 657     {
 658         /* Every molecule is a single charge group: no pbc required */
 659         pbcdxr = 0;
 660     }
 661
 662     if (cellSizeLimit > 0)
 663     {
 664         std::string maximumCells = "The maximum allowed number of cells is:";
 665         for (int d = 0; d < DIM; d++)
 666         {
 667             int nmax = static_cast<int>(ddbox.box_size[d] * ddbox.skew_fac[d] / cellSizeLimit);
 668             if (d >= ddbox.npbcdim && nmax < 2)
 669             {
 670                 nmax = 2;
 671             }
 672             if (d == ZZ && inhomogeneous_z(ir))
 673             {
 674                 nmax = 1;
 675             }
 676             maximumCells += gmx::formatString(" %c %d", 'X' + d, nmax);
 677         }
 678         GMX_LOG(mdlog.info).appendText(maximumCells);
 679     }
 680
 681     if (debug)
 682     {
 683         fprintf(debug, "Average nr of pbc_dx calls per atom %.2f\n", pbcdxr);
 684     }
 685
 686     /* Decompose numPPRanks in factors */
 687     std::vector<int> div;
 688     std::vector<int> mdiv;
 689     factorize(numPPRanks, &div, &mdiv);
 690
 691     gmx::IVec itry       = { 1, 1, 1 };
 692     gmx::IVec numDomains = { 0, 0, 0 };
 693     assign_factors(cellSizeLimit,
 694                    systemInfo.cutoff,
 695                    box,
 696                    ddbox,
 697                    mtop.natoms,
 698                    ir,
 699                    pbcdxr,
 700                    numRanksDoingPmeWork,
 701                    div.size(),
 702                    div.data(),
 703                    mdiv.data(),
 704                    &itry,
 705                    &numDomains);
 706
 707     return numDomains;
 708 }
 709
 710 real getDDGridSetupCellSizeLimit(const gmx::MDLogger& mdlog,
 711                                  const bool           bDynLoadBal,
 712                                  const real           dlb_scale,
 713                                  const t_inputrec&    ir,
 714                                  real                 systemInfoCellSizeLimit)
 715 {
 716     real cellSizeLimit = systemInfoCellSizeLimit;
 717
 718     /* Add a margin for DLB and/or pressure scaling */
 719     if (bDynLoadBal)
 720     {
 721         if (dlb_scale >= 1.0)
 722         {
 723             gmx_fatal(FARGS, "The value for option -dds should be smaller than 1");
 724         }
 725         GMX_LOG(mdlog.info)
 726                 .appendTextFormatted(
 727                         "Scaling the initial minimum size with 1/%g (option -dds) = %g", dlb_scale, 1 / dlb_scale);
 728         cellSizeLimit /= dlb_scale;
 729     }
 730     else if (ir.epc != PressureCoupling::No)
 731     {
 732         GMX_LOG(mdlog.info)
 733                 .appendTextFormatted(
 734                         "To account for pressure scaling, scaling the initial minimum size with %g",
 735                         DD_GRID_MARGIN_PRES_SCALE);
 736         cellSizeLimit *= DD_GRID_MARGIN_PRES_SCALE;
 737     }
 738
 739     return cellSizeLimit;
 740 }
 741 void checkForValidRankCountRequests(const int  numRanksRequested,
 742                                     const bool usingPme,
 743                                     const int  numPmeRanksRequested,
 744                                     const bool checkForLargePrimeFactors)
 745 {
 746     int numPPRanksRequested = numRanksRequested;
 747     if (usingPme && numPmeRanksRequested > 0)
 748     {
 749         numPPRanksRequested -= numPmeRanksRequested;
 750         if (numPmeRanksRequested > numPPRanksRequested)
 751         {
 752             gmx_fatal(FARGS,
 753                       "Cannot have %d separate PME ranks with only %d PP ranks, choose fewer or no "
 754                       "separate PME ranks",
 755                       numPmeRanksRequested,
 756                       numPPRanksRequested);
 757         }
 758     }
 759
 760     // Once the rank count is large enough, it becomes worth
 761     // suggesting improvements to the user.
 762     const int minPPRankCountToCheckForLargePrimeFactors = 13;
 763     if (checkForLargePrimeFactors && numPPRanksRequested >= minPPRankCountToCheckForLargePrimeFactors)
 764     {
 765         const int largestDivisor = largest_divisor(numPPRanksRequested);
 766         /* Check if the largest divisor is more than numPPRanks ^ (2/3) */
 767         if (largestDivisor * largestDivisor * largestDivisor > numPPRanksRequested * numPPRanksRequested)
 768         {
 769             gmx_fatal(FARGS,
 770                       "The number of ranks selected for particle-particle work (%d) "
 771                       "contains a large prime factor %d. In most cases this will lead to "
 772                       "bad performance. Choose a number with smaller prime factors or "
 773                       "set the decomposition (option -dd) manually.",
 774                       numPPRanksRequested,
 775                       largestDivisor);
 776         }
 777     }
 778 }
 779
 780 /*! \brief Return the number of PME-only ranks used by the simulation
 781  *
 782  * If the user did not choose a number, then decide for them. */
 783 static int getNumPmeOnlyRanksToUse(const gmx::MDLogger& mdlog,
 784                                    const DomdecOptions& options,
 785                                    const gmx_mtop_t&    mtop,
 786                                    const t_inputrec&    ir,
 787                                    const matrix         box,
 788                                    const int            numRanksRequested)
 789 {
 790     int         numPmeOnlyRanks = 0;
 791     const char* extraMessage    = "";
 792
 793     if (options.numCells[XX] > 0)
 794     {
 795         if (options.numPmeRanks >= 0)
 796         {
 797             // TODO mdrun should give a fatal error with a non-PME input file and -npme > 0
 798             numPmeOnlyRanks = options.numPmeRanks;
 799         }
 800         else
 801         {
 802             // When the DD grid is set explicitly and -npme is set to auto,
 803             // don't use PME ranks. We check later if the DD grid is
 804             // compatible with the total number of ranks.
 805             numPmeOnlyRanks = 0;
 806         }
 807     }
 808     else
 809     {
 810         if (!EEL_PME(ir.coulombtype))
 811         {
 812             // TODO mdrun should give a fatal error with a non-PME input file and -npme > 0
 813             numPmeOnlyRanks = 0;
 814         }
 815         else
 816         {
 817             if (options.numPmeRanks >= 0)
 818             {
 819                 numPmeOnlyRanks = options.numPmeRanks;
 820             }
 821             else
 822             {
 823                 // Controls the automated choice of when to use separate PME-only ranks.
 824                 const int minRankCountToDefaultToSeparatePmeRanks = 19;
 825
 826                 if (numRanksRequested < minRankCountToDefaultToSeparatePmeRanks)
 827                 {
 828                     numPmeOnlyRanks = 0;
 829                     extraMessage =
 830                             ", as there are too few total\n"
 831                             " ranks for efficient splitting";
 832                 }
 833                 else
 834                 {
 835                     numPmeOnlyRanks = guess_npme(mdlog, mtop, ir, box, numRanksRequested);
 836                     extraMessage    = ", as guessed by mdrun";
 837                 }
 838             }
 839         }
 840     }
 841     GMX_RELEASE_ASSERT(numPmeOnlyRanks <= numRanksRequested,
 842                        "Cannot have more PME ranks than total ranks");
 843     if (EEL_PME(ir.coulombtype))
 844     {
 845         GMX_LOG(mdlog.info).appendTextFormatted("Using %d separate PME ranks%s", numPmeOnlyRanks, extraMessage);
 846     }
 847
 848     return numPmeOnlyRanks;
 849 }
 850
 851 /*! \brief Sets the order of the DD dimensions, returns the number of DD dimensions */
 852 static int set_dd_dim(const gmx::IVec& numDDCells, const DDSettings& ddSettings, ivec* dims)
 853 {
 854     int ndim = 0;
 855     if (ddSettings.useDDOrderZYX)
 856     {
 857         /* Decomposition order z,y,x */
 858         for (int dim = DIM - 1; dim >= 0; dim--)
 859         {
 860             if (numDDCells[dim] > 1)
 861             {
 862                 (*dims)[ndim++] = dim;
 863             }
 864         }
 865     }
 866     else
 867     {
 868         /* Decomposition order x,y,z */
 869         for (int dim = 0; dim < DIM; dim++)
 870         {
 871             if (numDDCells[dim] > 1)
 872             {
 873                 (*dims)[ndim++] = dim;
 874             }
 875         }
 876     }
 877
 878     if (ndim == 0)
 879     {
 880         /* Set dim[0] to avoid extra checks on ndim in several places */
 881         (*dims)[0] = XX;
 882     }
 883
 884     return ndim;
 885 }
 886
 887 DDGridSetup getDDGridSetup(const gmx::MDLogger&           mdlog,
 888                            DDRole                         ddRole,
 889                            MPI_Comm                       communicator,
 890                            const int                      numRanksRequested,
 891                            const DomdecOptions&           options,
 892                            const DDSettings&              ddSettings,
 893                            const DDSystemInfo&            systemInfo,
 894                            const real                     cellSizeLimit,
 895                            const gmx_mtop_t&              mtop,
 896                            const t_inputrec&              ir,
 897                            const matrix                   box,
 898                            gmx::ArrayRef<const gmx::RVec> xGlobal,
 899                            gmx_ddbox_t*                   ddbox)
 900 {
 901     int numPmeOnlyRanks = getNumPmeOnlyRanksToUse(mdlog, options, mtop, ir, box, numRanksRequested);
 902
 903     gmx::IVec numDomains;
 904     if (options.numCells[XX] > 0)
 905     {
 906         numDomains                      = gmx::IVec(options.numCells);
 907         const ivec numDomainsLegacyIvec = { numDomains[XX], numDomains[YY], numDomains[ZZ] };
 908         set_ddbox_cr(ddRole, communicator, &numDomainsLegacyIvec, ir, box, xGlobal, ddbox);
 909     }
 910     else
 911     {
 912         set_ddbox_cr(ddRole, communicator, nullptr, ir, box, xGlobal, ddbox);
 913
 914         if (ddRole == DDRole::Master)
 915         {
 916             numDomains = optimizeDDCells(
 917                     mdlog, numRanksRequested, numPmeOnlyRanks, cellSizeLimit, mtop, box, *ddbox, ir, systemInfo);
 918         }
 919     }
 920
 921     /* Communicate the information set by the master to all ranks */
 922     gmx_bcast(sizeof(numDomains), numDomains, communicator);
 923     if (EEL_PME(ir.coulombtype))
 924     {
 925         gmx_bcast(sizeof(numPmeOnlyRanks), &numPmeOnlyRanks, communicator);
 926     }
 927
 928     DDGridSetup ddGridSetup;
 929     ddGridSetup.numPmeOnlyRanks = numPmeOnlyRanks;
 930     ddGridSetup.numDomains[XX]  = numDomains[XX];
 931     ddGridSetup.numDomains[YY]  = numDomains[YY];
 932     ddGridSetup.numDomains[ZZ]  = numDomains[ZZ];
 933     ddGridSetup.numDDDimensions = set_dd_dim(numDomains, ddSettings, &ddGridSetup.ddDimensions);
 934
 935     return ddGridSetup;
 936 }