src/gromacs/mdlib/domdec_setup.c

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