src/gromacs/mdlib/domdec_setup.c

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