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