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