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