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