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