src/mdlib/domdec_setup.c

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