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