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