src/gromacs/domdec/domdec_setup.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2008,2009,2010,2011,2012,2013,2014,2015,2017,2018,2019, 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 /*! \internal \file
  37  *
  38  * \brief This file defines functions used by the domdec module
  39  * in its initial setup phase.
  40  *
  41  * \author Berk Hess <hess@kth.se>
  42  * \ingroup module_domdec
  43  */
  44
  45 #include "gmxpre.h"
  46
  47 #include "domdec_setup.h"
  48
  49 #include <cassert>
  50 #include <cinttypes>
  51 #include <cmath>
  52 #include <cstdio>
  53
  54 #include "gromacs/domdec/domdec.h"
  55 #include "gromacs/domdec/domdec_struct.h"
  56 #include "gromacs/domdec/options.h"
  57 #include "gromacs/ewald/pme.h"
  58 #include "gromacs/gmxlib/network.h"
  59 #include "gromacs/math/utilities.h"
  60 #include "gromacs/math/vec.h"
  61 #include "gromacs/mdlib/perf_est.h"
  62 #include "gromacs/mdtypes/commrec.h"
  63 #include "gromacs/mdtypes/inputrec.h"
  64 #include "gromacs/mdtypes/md_enums.h"
  65 #include "gromacs/pbcutil/pbc.h"
  66 #include "gromacs/topology/topology.h"
  67 #include "gromacs/utility/fatalerror.h"
  68 #include "gromacs/utility/logger.h"
  69 #include "gromacs/utility/stringutil.h"
  70
  71 #include "box.h"
  72 #include "domdec_internal.h"
  73 #include "utility.h"
  74
  75 // TODO remove this when moving domdec into gmx namespace
  76 using gmx::DomdecOptions;
  77
  78 /*! \brief Margin for setting up the DD grid */
  79 #define DD_GRID_MARGIN_PRES_SCALE 1.05
  80
  81 /*! \brief Factorize \p n.
  82  *
  83  * \param[in]    n     Value to factorize
  84  * \param[out]   fac   Vector of factors (to be allocated in this function)
  85  * \param[out]   mfac  Vector with the number of times each factor repeats in the factorization (to be allocated in this function)
  86  */
  87 static void factorize(int               n,
  88                       std::vector<int> *fac,
  89                       std::vector<int> *mfac)
  90 {
  91     if (n <= 0)
  92     {
  93         gmx_fatal(FARGS, "Can only factorize positive integers.");
  94     }
  95
  96     /* Decompose n in factors */
  97     fac->clear();
  98     mfac->clear();
  99     int d = 2;
 100     while (n > 1)
 101     {
 102         while (n % d == 0)
 103         {
 104             if (fac->empty() || fac->back() != d)
 105             {
 106                 fac->push_back(d);
 107                 mfac->push_back(1);
 108             }
 109             else
 110             {
 111                 mfac->back()++;
 112             }
 113             n /= d;
 114         }
 115         d++;
 116     }
 117 }
 118
 119 /*! \brief Find largest divisor of \p n smaller than \p n*/
 120 static int largest_divisor(int n)
 121 {
 122     std::vector<int> div;
 123     std::vector<int> mdiv;
 124     factorize(n, &div, &mdiv);
 125
 126     return div.back();
 127 }
 128
 129 /*! \brief Compute largest common divisor of \p n1 and \b n2 */
 130 static int lcd(int n1, int n2)
 131 {
 132     int d, i;
 133
 134     d = 1;
 135     for (i = 2; (i <= n1 && i <= n2); i++)
 136     {
 137         if (n1 % i == 0 && n2 % i == 0)
 138         {
 139             d = i;
 140         }
 141     }
 142
 143     return d;
 144 }
 145
 146 /*! \brief Returns TRUE when there are enough PME ranks for the ratio */
 147 static gmx_bool fits_pme_ratio(int nrank_tot, int nrank_pme, float ratio)
 148 {
 149     return (static_cast<double>(nrank_pme)/static_cast<double>(nrank_tot) > 0.95*ratio);
 150 }
 151
 152 /*! \brief Returns TRUE when npme out of ntot ranks doing PME is expected to give reasonable performance */
 153 static gmx_bool fits_pp_pme_perf(int ntot, int npme, float ratio)
 154 {
 155     std::vector<int> div;
 156     std::vector<int> mdiv;
 157     factorize(ntot - npme, &div, &mdiv);
 158
 159     int npp_root3  = gmx::roundToInt(std::cbrt(ntot - npme));
 160     int npme_root2 = gmx::roundToInt(std::sqrt(static_cast<double>(npme)));
 161
 162     /* The check below gives a reasonable division:
 163      * factor 5 allowed at 5 or more PP ranks,
 164      * factor 7 allowed at 49 or more PP ranks.
 165      */
 166     if (div.back() > 3 + npp_root3)
 167     {
 168         return FALSE;
 169     }
 170
 171     /* Check if the number of PP and PME ranks have a reasonable sized
 172      * denominator in common, such that we can use 2D PME decomposition
 173      * when required (which requires nx_pp == nx_pme).
 174      * The factor of 2 allows for a maximum ratio of 2^2=4
 175      * between nx_pme and ny_pme.
 176      */
 177     if (lcd(ntot - npme, npme)*2 < npme_root2)
 178     {
 179         return FALSE;
 180     }
 181
 182     /* Does this division gives a reasonable PME load? */
 183     return fits_pme_ratio(ntot, npme, ratio);
 184 }
 185
 186 /*! \brief Make a guess for the number of PME ranks to use. */
 187 static int guess_npme(const gmx::MDLogger &mdlog,
 188                       const gmx_mtop_t    &mtop,
 189                       const t_inputrec    &ir,
 190                       const matrix         box,
 191                       int                  nrank_tot)
 192 {
 193     float      ratio;
 194     int        npme;
 195
 196     ratio = pme_load_estimate(mtop, ir, box);
 197
 198     GMX_LOG(mdlog.info).appendTextFormatted(
 199             "Guess for relative PME load: %.2f", ratio);
 200
 201     /* We assume the optimal rank ratio is close to the load ratio.
 202      * The communication load is neglected,
 203      * but (hopefully) this will balance out between PP and PME.
 204      */
 205
 206     if (!fits_pme_ratio(nrank_tot, nrank_tot/2, ratio))
 207     {
 208         /* We would need more than nrank_tot/2 PME only nodes,
 209          * which is not possible. Since the PME load is very high,
 210          * we will not loose much performance when all ranks do PME.
 211          */
 212
 213         return 0;
 214     }
 215
 216     /* First try to find npme as a factor of nrank_tot up to nrank_tot/3.
 217      * We start with a minimum PME node fraction of 1/16
 218      * and avoid ratios which lead to large prime factors in nnodes-npme.
 219      */
 220     npme = (nrank_tot + 15)/16;
 221     while (npme <= nrank_tot/3)
 222     {
 223         if (nrank_tot % npme == 0)
 224         {
 225             /* Note that fits_perf might change the PME grid,
 226              * in the current implementation it does not.
 227              */
 228             if (fits_pp_pme_perf(nrank_tot, npme, ratio))
 229             {
 230                 break;
 231             }
 232         }
 233         npme++;
 234     }
 235     if (npme > nrank_tot/3)
 236     {
 237         /* Try any possible number for npme */
 238         npme = 1;
 239         while (npme <= nrank_tot/2)
 240         {
 241             /* Note that fits_perf may change the PME grid */
 242             if (fits_pp_pme_perf(nrank_tot, npme, ratio))
 243             {
 244                 break;
 245             }
 246             npme++;
 247         }
 248     }
 249     if (npme > nrank_tot/2)
 250     {
 251         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"
 252                   "Use the -npme option of mdrun or change the number of ranks or the PME grid dimensions, see the manual for details.",
 253                   ratio, gmx::roundToInt(0.95*ratio*nrank_tot), nrank_tot/2, ir.nkx, ir.nky);
 254     }
 255     else
 256     {
 257         GMX_LOG(mdlog.info).appendTextFormatted(
 258                 "Will use %d particle-particle and %d PME only ranks\n"
 259                 "This is a guess, check the performance at the end of the log file",
 260                 nrank_tot - npme, npme);
 261     }
 262
 263     return npme;
 264 }
 265
 266 /*! \brief Return \p n divided by \p f rounded up to the next integer. */
 267 static int div_up(int n, int f)
 268 {
 269     return (n + f - 1)/f;
 270 }
 271
 272 real comm_box_frac(const gmx::IVec &dd_nc, real cutoff, const gmx_ddbox_t &ddbox)
 273 {
 274     int  i, j, k;
 275     rvec nw;
 276     real comm_vol;
 277
 278     for (i = 0; i < DIM; i++)
 279     {
 280         real bt = ddbox.box_size[i]*ddbox.skew_fac[i];
 281         nw[i] = dd_nc[i]*cutoff/bt;
 282     }
 283
 284     comm_vol = 0;
 285     for (i = 0; i < DIM; i++)
 286     {
 287         if (dd_nc[i] > 1)
 288         {
 289             comm_vol += nw[i];
 290             for (j = i+1; j < DIM; j++)
 291             {
 292                 if (dd_nc[j] > 1)
 293                 {
 294                     comm_vol += nw[i]*nw[j]*M_PI/4;
 295                     for (k = j+1; k < DIM; k++)
 296                     {
 297                         if (dd_nc[k] > 1)
 298                         {
 299                             comm_vol += nw[i]*nw[j]*nw[k]*M_PI/6;
 300                         }
 301                     }
 302                 }
 303             }
 304         }
 305     }
 306
 307     return comm_vol;
 308 }
 309
 310 /*! \brief Return whether the DD inhomogeneous in the z direction */
 311 static gmx_bool inhomogeneous_z(const t_inputrec &ir)
 312 {
 313     return ((EEL_PME(ir.coulombtype) || ir.coulombtype == eelEWALD) &&
 314             ir.ePBC == epbcXYZ && ir.ewald_geometry == eewg3DC);
 315 }
 316
 317 /*! \brief Estimate cost of PME FFT communication
 318  *
 319  * This only takes the communication into account and not imbalance
 320  * in the calculation. But the imbalance in communication and calculation
 321  * are similar and therefore these formulas also prefer load balance
 322  * in the FFT and pme_solve calculation.
 323  */
 324 static float comm_pme_cost_vol(int npme, int a, int b, int c)
 325 {
 326     /* We use a float here, since an integer might overflow */
 327     float comm_vol;
 328
 329     comm_vol  = npme - 1;
 330     comm_vol *= npme;
 331     comm_vol *= div_up(a, npme);
 332     comm_vol *= div_up(b, npme);
 333     comm_vol *= c;
 334
 335     return comm_vol;
 336 }
 337
 338 /*! \brief Estimate cost of communication for a possible domain decomposition. */
 339 static float comm_cost_est(real limit, real cutoff,
 340                            const matrix box, const gmx_ddbox_t &ddbox,
 341                            int natoms, const t_inputrec &ir,
 342                            float pbcdxr,
 343                            int npme_tot, const gmx::IVec &nc)
 344 {
 345     gmx::IVec npme = {1, 1, 1};
 346     int       i, j, nk, overlap;
 347     rvec      bt;
 348     float     comm_vol, comm_vol_xf, comm_pme, cost_pbcdx;
 349     /* This is the cost of a pbc_dx call relative to the cost
 350      * of communicating the coordinate and force of an atom.
 351      * This will be machine dependent.
 352      * These factors are for x86 with SMP or Infiniband.
 353      */
 354     float pbcdx_rect_fac = 0.1;
 355     float pbcdx_tric_fac = 0.2;
 356     float temp;
 357
 358     /* Check the DD algorithm restrictions */
 359     if ((ir.ePBC == epbcXY && ir.nwall < 2 && nc[ZZ] > 1) ||
 360         (ir.ePBC == epbcSCREW && (nc[XX] == 1 || nc[YY] > 1 || nc[ZZ] > 1)))
 361     {
 362         return -1;
 363     }
 364
 365     if (inhomogeneous_z(ir) && nc[ZZ] > 1)
 366     {
 367         return -1;
 368     }
 369
 370     assert(ddbox.npbcdim <= DIM);
 371
 372     /* Check if the triclinic requirements are met */
 373     for (i = 0; i < DIM; i++)
 374     {
 375         for (j = i+1; j < ddbox.npbcdim; j++)
 376         {
 377             if (box[j][i] != 0 || ir.deform[j][i] != 0 ||
 378                 (ir.epc != epcNO && ir.compress[j][i] != 0))
 379             {
 380                 if (nc[j] > 1 && nc[i] == 1)
 381                 {
 382                     return -1;
 383                 }
 384             }
 385         }
 386     }
 387
 388     for (i = 0; i < DIM; i++)
 389     {
 390         bt[i] = ddbox.box_size[i]*ddbox.skew_fac[i];
 391
 392         /* Without PBC and with 2 cells, there are no lower limits on the cell size */
 393         if (!(i >= ddbox.npbcdim && nc[i] <= 2) && bt[i] < nc[i]*limit)
 394         {
 395             return -1;
 396         }
 397         /* With PBC, check if the cut-off fits in nc[i]-1 cells */
 398         if (i < ddbox.npbcdim && nc[i] > 1 && (nc[i] - 1)*bt[i] < nc[i]*cutoff)
 399         {
 400             return -1;
 401         }
 402     }
 403
 404     if (npme_tot > 1)
 405     {
 406         /* The following choices should match those
 407          * in init_domain_decomposition in domdec.c.
 408          */
 409         if (nc[XX] == 1 && nc[YY] > 1)
 410         {
 411             npme[XX] = 1;
 412             npme[YY] = npme_tot;
 413         }
 414         else if (nc[YY] == 1)
 415         {
 416             npme[XX] = npme_tot;
 417             npme[YY] = 1;
 418         }
 419         else
 420         {
 421             /* Will we use 1D or 2D PME decomposition? */
 422             npme[XX] = (npme_tot % nc[XX] == 0) ? nc[XX] : npme_tot;
 423             npme[YY] = npme_tot/npme[XX];
 424         }
 425     }
 426
 427     if (EEL_PME(ir.coulombtype) || EVDW_PME(ir.vdwtype))
 428     {
 429         /* Check the PME grid restrictions.
 430          * Currently these can only be invalid here with too few grid lines
 431          * along the x dimension per rank doing PME.
 432          */
 433         int npme_x = (npme_tot > 1 ? npme[XX] : nc[XX]);
 434
 435         /* Currently we don't have the OpenMP thread count available here.
 436          * But with threads we have only tighter restrictions and it's
 437          * probably better anyhow to avoid settings where we need to reduce
 438          * grid lines over multiple ranks, as the thread check will do.
 439          */
 440         bool useThreads     = true;
 441         bool errorsAreFatal = false;
 442         if (!gmx_pme_check_restrictions(ir.pme_order, ir.nkx, ir.nky, ir.nkz,
 443                                         npme_x, useThreads, errorsAreFatal))
 444         {
 445             return -1;
 446         }
 447     }
 448
 449     /* When two dimensions are (nearly) equal, use more cells
 450      * for the smallest index, so the decomposition does not
 451      * depend sensitively on the rounding of the box elements.
 452      */
 453     for (i = 0; i < DIM; i++)
 454     {
 455         for (j = i+1; j < DIM; j++)
 456         {
 457             /* Check if the box size is nearly identical,
 458              * in that case we prefer nx > ny  and ny > nz.
 459              */
 460             if (std::fabs(bt[j] - bt[i]) < 0.01*bt[i] && nc[j] > nc[i])
 461             {
 462                 /* The XX/YY check is a bit compact. If nc[YY]==npme[YY]
 463                  * this means the swapped nc has nc[XX]==npme[XX],
 464                  * and we can also swap X and Y for PME.
 465                  */
 466                 /* Check if dimension i and j are equivalent for PME.
 467                  * For x/y: if nc[YY]!=npme[YY], we can not swap x/y
 468                  * For y/z: we can not have PME decomposition in z
 469                  */
 470                 if (npme_tot <= 1 ||
 471                     !((i == XX && j == YY && nc[YY] != npme[YY]) ||
 472                       (i == YY && j == ZZ && npme[YY] > 1)))
 473                 {
 474                     return -1;
 475                 }
 476             }
 477         }
 478     }
 479
 480     /* This function determines only half of the communication cost.
 481      * All PP, PME and PP-PME communication is symmetric
 482      * and the "back"-communication cost is identical to the forward cost.
 483      */
 484
 485     comm_vol = comm_box_frac(nc, cutoff, ddbox);
 486
 487     comm_pme = 0;
 488     for (i = 0; i < 2; i++)
 489     {
 490         /* Determine the largest volume for PME x/f redistribution */
 491         if (nc[i] % npme[i] != 0)
 492         {
 493             if (nc[i] > npme[i])
 494             {
 495                 comm_vol_xf = (npme[i] == 2 ? 1.0/3.0 : 0.5);
 496             }
 497             else
 498             {
 499                 comm_vol_xf = 1.0 - lcd(nc[i], npme[i])/static_cast<double>(npme[i]);
 500             }
 501             comm_pme += 3*natoms*comm_vol_xf;
 502         }
 503
 504         /* Grid overlap communication */
 505         if (npme[i] > 1)
 506         {
 507             nk        = (i == 0 ? ir.nkx : ir.nky);
 508             overlap   = (nk % npme[i] == 0 ? ir.pme_order-1 : ir.pme_order);
 509             temp      = npme[i];
 510             temp     *= overlap;
 511             temp     *= ir.nkx;
 512             temp     *= ir.nky;
 513             temp     *= ir.nkz;
 514             temp     /= nk;
 515             comm_pme += temp;
 516 /* Old line comm_pme += npme[i]*overlap*ir.nkx*ir.nky*ir.nkz/nk; */
 517         }
 518     }
 519
 520     comm_pme += comm_pme_cost_vol(npme[YY], ir.nky, ir.nkz, ir.nkx);
 521     comm_pme += comm_pme_cost_vol(npme[XX], ir.nkx, ir.nky, ir.nkz);
 522
 523     /* Add cost of pbc_dx for bondeds */
 524     cost_pbcdx = 0;
 525     if ((nc[XX] == 1 || nc[YY] == 1) || (nc[ZZ] == 1 && ir.ePBC != epbcXY))
 526     {
 527         if ((ddbox.tric_dir[XX] && nc[XX] == 1) ||
 528             (ddbox.tric_dir[YY] && nc[YY] == 1))
 529         {
 530             cost_pbcdx = pbcdxr*pbcdx_tric_fac;
 531         }
 532         else
 533         {
 534             cost_pbcdx = pbcdxr*pbcdx_rect_fac;
 535         }
 536     }
 537
 538     if (debug)
 539     {
 540         fprintf(debug,
 541                 "nc %2d %2d %2d %2d %2d vol pp %6.4f pbcdx %6.4f pme %9.3e tot %9.3e\n",
 542                 nc[XX], nc[YY], nc[ZZ], npme[XX], npme[YY],
 543                 comm_vol, cost_pbcdx, comm_pme/(3*natoms),
 544                 comm_vol + cost_pbcdx + comm_pme/(3*natoms));
 545     }
 546
 547     return 3*natoms*(comm_vol + cost_pbcdx) + comm_pme;
 548 }
 549
 550 /*! \brief Assign penalty factors to possible domain decompositions,
 551  * based on the estimated communication costs. */
 552 static void assign_factors(const real limit, const bool request1D,
 553                            const real cutoff,
 554                            const matrix box, const gmx_ddbox_t &ddbox,
 555                            int natoms, const t_inputrec &ir,
 556                            float pbcdxr, int npme,
 557                            int ndiv, const int *div, const int *mdiv,
 558                            gmx::IVec *irTryPtr,
 559                            gmx::IVec *opt)
 560 {
 561     int        x, y, i;
 562     float      ce;
 563     gmx::IVec &ir_try = *irTryPtr;
 564
 565     if (ndiv == 0)
 566     {
 567         const int  maxDimensionSize             = std::max(ir_try[XX], std::max(ir_try[YY], ir_try[ZZ]));
 568         const int  productOfDimensionSizes      = ir_try[XX]*ir_try[YY]*ir_try[ZZ];
 569         const bool decompositionHasOneDimension = (maxDimensionSize == productOfDimensionSizes);
 570         if (request1D && !decompositionHasOneDimension)
 571         {
 572             /* We requested 1D DD, but got multiple dimensions */
 573             return;
 574         }
 575
 576         ce = comm_cost_est(limit, cutoff, box, ddbox,
 577                            natoms, ir, pbcdxr, npme, ir_try);
 578         if (ce >= 0 && ((*opt)[XX] == 0 ||
 579                         ce < comm_cost_est(limit, cutoff, box, ddbox,
 580                                            natoms, ir, pbcdxr,
 581                                            npme, *opt)))
 582         {
 583             *opt = ir_try;
 584         }
 585
 586         return;
 587     }
 588
 589     for (x = mdiv[0]; x >= 0; x--)
 590     {
 591         for (i = 0; i < x; i++)
 592         {
 593             ir_try[XX] *= div[0];
 594         }
 595         for (y = mdiv[0]-x; y >= 0; y--)
 596         {
 597             for (i = 0; i < y; i++)
 598             {
 599                 ir_try[YY] *= div[0];
 600             }
 601             for (i = 0; i < mdiv[0]-x-y; i++)
 602             {
 603                 ir_try[ZZ] *= div[0];
 604             }
 605
 606             /* recurse */
 607             assign_factors(limit, request1D,
 608                            cutoff, box, ddbox, natoms, ir, pbcdxr, npme,
 609                            ndiv-1, div+1, mdiv+1, irTryPtr, opt);
 610
 611             for (i = 0; i < mdiv[0]-x-y; i++)
 612             {
 613                 ir_try[ZZ] /= div[0];
 614             }
 615             for (i = 0; i < y; i++)
 616             {
 617                 ir_try[YY] /= div[0];
 618             }
 619         }
 620         for (i = 0; i < x; i++)
 621         {
 622             ir_try[XX] /= div[0];
 623         }
 624     }
 625 }
 626
 627 /*! \brief Determine the optimal distribution of DD cells for the
 628  * simulation system and number of MPI ranks
 629  *
 630  * \returns The optimal grid cell choice. The latter will contain all
 631  *          zeros if no valid cell choice exists. */
 632 static gmx::IVec
 633 optimizeDDCells(const gmx::MDLogger &mdlog,
 634                 const int            numRanksRequested,
 635                 const int            numPmeOnlyRanks,
 636                 const real           cellSizeLimit,
 637                 const bool           request1DAnd1Pulse,
 638                 const gmx_mtop_t    &mtop,
 639                 const matrix         box,
 640                 const gmx_ddbox_t   &ddbox,
 641                 const t_inputrec    &ir,
 642                 const DDSystemInfo  &systemInfo)
 643 {
 644     double    pbcdxr;
 645
 646     const int numPPRanks = numRanksRequested - numPmeOnlyRanks;
 647
 648     GMX_LOG(mdlog.info).appendTextFormatted(
 649             "Optimizing the DD grid for %d cells with a minimum initial size of %.3f nm",
 650             numPPRanks, cellSizeLimit);
 651     if (inhomogeneous_z(ir))
 652     {
 653         GMX_LOG(mdlog.info).appendTextFormatted(
 654                 "Ewald_geometry=%s: assuming inhomogeneous particle distribution in z, will not decompose in z.",
 655                 eewg_names[ir.ewald_geometry]);
 656     }
 657
 658
 659     // For cost estimates, we need the number of ranks doing PME work,
 660     // which is the number of PP ranks when not using separate
 661     // PME-only ranks.
 662     const int numRanksDoingPmeWork = (EEL_PME(ir.coulombtype) ?
 663                                       ((numPmeOnlyRanks > 0) ? numPmeOnlyRanks : numPPRanks) :
 664                                       0);
 665
 666     if (systemInfo.haveInterDomainBondeds)
 667     {
 668         /* If we can skip PBC for distance calculations in plain-C bondeds,
 669          * we can save some time (e.g. 3D DD with pbc=xyz).
 670          * Here we ignore SIMD bondeds as they always do (fast) PBC.
 671          */
 672         count_bonded_distances(mtop, ir, &pbcdxr, nullptr);
 673         pbcdxr /= static_cast<double>(mtop.natoms);
 674     }
 675     else
 676     {
 677         /* Every molecule is a single charge group: no pbc required */
 678         pbcdxr = 0;
 679     }
 680
 681     if (cellSizeLimit > 0)
 682     {
 683         std::string maximumCells = "The maximum allowed number of cells is:";
 684         for (int d = 0; d < DIM; d++)
 685         {
 686             int nmax = static_cast<int>(ddbox.box_size[d]*ddbox.skew_fac[d]/cellSizeLimit);
 687             if (d >= ddbox.npbcdim && nmax < 2)
 688             {
 689                 nmax = 2;
 690             }
 691             if (d == ZZ && inhomogeneous_z(ir))
 692             {
 693                 nmax = 1;
 694             }
 695             maximumCells += gmx::formatString(" %c %d", 'X' + d, nmax);
 696         }
 697         GMX_LOG(mdlog.info).appendText(maximumCells);
 698     }
 699
 700     if (debug)
 701     {
 702         fprintf(debug, "Average nr of pbc_dx calls per atom %.2f\n", pbcdxr);
 703     }
 704
 705     /* Decompose numPPRanks in factors */
 706     std::vector<int> div;
 707     std::vector<int> mdiv;
 708     factorize(numPPRanks, &div, &mdiv);
 709
 710     gmx::IVec itry       = { 1, 1, 1 };
 711     gmx::IVec numDomains = { 0, 0, 0 };
 712     assign_factors(cellSizeLimit, request1DAnd1Pulse,
 713                    systemInfo.cutoff, box, ddbox, mtop.natoms, ir, pbcdxr,
 714                    numRanksDoingPmeWork, div.size(), div.data(), mdiv.data(), &itry, &numDomains);
 715
 716     return numDomains;
 717 }
 718
 719 real
 720 getDDGridSetupCellSizeLimit(const gmx::MDLogger &mdlog,
 721                             const bool           request1DAnd1Pulse,
 722                             const bool           bDynLoadBal,
 723                             const real           dlb_scale,
 724                             const t_inputrec    &ir,
 725                             const real           systemInfoCellSizeLimit)
 726 {
 727     real cellSizeLimit = systemInfoCellSizeLimit;
 728     if (request1DAnd1Pulse)
 729     {
 730         cellSizeLimit = std::max(cellSizeLimit, ir.rlist);
 731     }
 732
 733     /* Add a margin for DLB and/or pressure scaling */
 734     if (bDynLoadBal)
 735     {
 736         if (dlb_scale >= 1.0)
 737         {
 738             gmx_fatal(FARGS, "The value for option -dds should be smaller than 1");
 739         }
 740         GMX_LOG(mdlog.info).appendTextFormatted(
 741                 "Scaling the initial minimum size with 1/%g (option -dds) = %g",
 742                 dlb_scale, 1/dlb_scale);
 743         cellSizeLimit /= dlb_scale;
 744     }
 745     else if (ir.epc != epcNO)
 746     {
 747         GMX_LOG(mdlog.info).appendTextFormatted(
 748                 "To account for pressure scaling, scaling the initial minimum size with %g",
 749                 DD_GRID_MARGIN_PRES_SCALE);
 750         cellSizeLimit *= DD_GRID_MARGIN_PRES_SCALE;
 751     }
 752
 753     return cellSizeLimit;
 754 }
 755 void
 756 checkForValidRankCountRequests(const int  numRanksRequested,
 757                                const bool usingPme,
 758                                const int  numPmeRanksRequested)
 759 {
 760     int numPPRanksRequested = numRanksRequested;
 761     if (usingPme)
 762     {
 763         numPPRanksRequested -= numPmeRanksRequested;
 764         if (numPmeRanksRequested > numPPRanksRequested)
 765         {
 766             gmx_fatal(FARGS,
 767                       "Cannot have %d separate PME ranks with only %d PP ranks, choose fewer or no separate PME ranks",
 768                       numPmeRanksRequested, numPPRanksRequested);
 769         }
 770     }
 771
 772     // Once the rank count is large enough, it becomes worth
 773     // suggesting improvements to the user.
 774     const int minPPRankCountToCheckForLargePrimeFactors = 13;
 775     if (numPPRanksRequested >= minPPRankCountToCheckForLargePrimeFactors)
 776     {
 777         const int largestDivisor = largest_divisor(numPPRanksRequested);
 778         /* Check if the largest divisor is more than numPPRanks ^ (2/3) */
 779         if (largestDivisor*largestDivisor*largestDivisor >
 780             numPPRanksRequested*numPPRanksRequested)
 781         {
 782             gmx_fatal(FARGS, "The number of ranks selected for particle-particle work (%d) "
 783                       "contains a large prime factor %d. In most cases this will lead to "
 784                       "bad performance. Choose a number with smaller prime factors or "
 785                       "set the decomposition (option -dd) manually.",
 786                       numPPRanksRequested, largestDivisor);
 787         }
 788     }
 789 }
 790
 791 /*! \brief Return the number of PME-only ranks used by the simulation
 792  *
 793  * If the user did not choose a number, then decide for them. */
 794 static int
 795 getNumPmeOnlyRanksToUse(const gmx::MDLogger &mdlog,
 796                         const DomdecOptions &options,
 797                         const gmx_mtop_t    &mtop,
 798                         const t_inputrec    &ir,
 799                         const matrix         box,
 800                         const int            numRanksRequested)
 801 {
 802     int         numPmeOnlyRanks;
 803     const char *extraMessage = "";
 804
 805     if (options.numCells[XX] > 0)
 806     {
 807         if (options.numPmeRanks >= 0)
 808         {
 809             // TODO mdrun should give a fatal error with a non-PME input file and -npme > 0
 810             numPmeOnlyRanks = options.numPmeRanks;
 811         }
 812         else
 813         {
 814             // When the DD grid is set explicitly and -npme is set to auto,
 815             // don't use PME ranks. We check later if the DD grid is
 816             // compatible with the total number of ranks.
 817             numPmeOnlyRanks = 0;
 818         }
 819     }
 820     else
 821     {
 822         if (!EEL_PME(ir.coulombtype))
 823         {
 824             // TODO mdrun should give a fatal error with a non-PME input file and -npme > 0
 825             numPmeOnlyRanks = 0;
 826         }
 827         else
 828         {
 829             if (options.numPmeRanks >= 0)
 830             {
 831                 numPmeOnlyRanks = options.numPmeRanks;
 832             }
 833             else
 834             {
 835                 // Controls the automated choice of when to use separate PME-only ranks.
 836                 const int minRankCountToDefaultToSeparatePmeRanks = 19;
 837
 838                 if (numRanksRequested < minRankCountToDefaultToSeparatePmeRanks)
 839                 {
 840                     numPmeOnlyRanks = 0;
 841                     extraMessage    = ", as there are too few total\n"
 842                         " ranks for efficient splitting";
 843                 }
 844                 else
 845                 {
 846                     numPmeOnlyRanks = guess_npme(mdlog, mtop, ir, box, numRanksRequested);
 847                     extraMessage    = ", as guessed by mdrun";
 848                 }
 849             }
 850         }
 851
 852     }
 853     GMX_RELEASE_ASSERT(numPmeOnlyRanks <= numRanksRequested,
 854                        "Cannot have more PME ranks than total ranks");
 855     if (EEL_PME(ir.coulombtype))
 856     {
 857         GMX_LOG(mdlog.info).appendTextFormatted
 858             ("Using %d separate PME ranks%s", numPmeOnlyRanks, extraMessage);
 859     }
 860
 861     return numPmeOnlyRanks;
 862 }
 863
 864 /*! \brief Sets the order of the DD dimensions, returns the number of DD dimensions */
 865 static int set_dd_dim(const gmx::IVec  &numDDCells,
 866                       const DDSettings &ddSettings,
 867                       ivec             *dims)
 868 {
 869     int ndim = 0;
 870     if (ddSettings.useDDOrderZYX)
 871     {
 872         /* Decomposition order z,y,x */
 873         for (int dim = DIM - 1; dim >= 0; dim--)
 874         {
 875             if (numDDCells[dim] > 1)
 876             {
 877                 (*dims)[ndim++] = dim;
 878             }
 879         }
 880     }
 881     else
 882     {
 883         /* Decomposition order x,y,z */
 884         for (int dim = 0; dim < DIM; dim++)
 885         {
 886             if (numDDCells[dim] > 1)
 887             {
 888                 (*dims)[ndim++] = dim;
 889             }
 890         }
 891     }
 892
 893     if (ndim == 0)
 894     {
 895         /* Set dim[0] to avoid extra checks on ndim in several places */
 896         (*dims)[0] = XX;
 897     }
 898
 899     return ndim;
 900 }
 901
 902 DDGridSetup
 903 getDDGridSetup(const gmx::MDLogger           &mdlog,
 904                const t_commrec               *cr,
 905                const int                      numRanksRequested,
 906                const DomdecOptions           &options,
 907                const DDSettings              &ddSettings,
 908                const DDSystemInfo            &systemInfo,
 909                const real                     cellSizeLimit,
 910                const gmx_mtop_t              &mtop,
 911                const t_inputrec              &ir,
 912                const matrix                   box,
 913                gmx::ArrayRef<const gmx::RVec> xGlobal,
 914                gmx_ddbox_t                   *ddbox)
 915 {
 916     int       numPmeOnlyRanks = getNumPmeOnlyRanksToUse(mdlog, options, mtop, ir, box, numRanksRequested);
 917
 918     gmx::IVec numDomains;
 919     if (options.numCells[XX] > 0)
 920     {
 921         numDomains = gmx::IVec(options.numCells);
 922         const ivec numDomainsLegacyIvec = { numDomains[XX], numDomains[YY], numDomains[ZZ] };
 923         set_ddbox_cr(*cr, &numDomainsLegacyIvec, ir, box, xGlobal, ddbox);
 924     }
 925     else
 926     {
 927         set_ddbox_cr(*cr, nullptr, ir, box, xGlobal, ddbox);
 928
 929         if (MASTER(cr))
 930         {
 931             numDomains =
 932                 optimizeDDCells(mdlog, numRanksRequested, numPmeOnlyRanks,
 933                                 cellSizeLimit,
 934                                 ddSettings.request1DAnd1Pulse,
 935                                 mtop, box, *ddbox, ir,
 936                                 systemInfo);
 937         }
 938     }
 939
 940     /* Communicate the information set by the master to all ranks */
 941     gmx_bcast(sizeof(numDomains), numDomains, cr);
 942     if (EEL_PME(ir.coulombtype))
 943     {
 944         gmx_bcast(sizeof(numPmeOnlyRanks), &numPmeOnlyRanks, cr);
 945     }
 946
 947     DDGridSetup ddGridSetup;
 948     ddGridSetup.numPmeOnlyRanks = numPmeOnlyRanks;
 949     ddGridSetup.numDomains[XX]  = numDomains[XX];
 950     ddGridSetup.numDomains[YY]  = numDomains[YY];
 951     ddGridSetup.numDomains[ZZ]  = numDomains[ZZ];
 952     ddGridSetup.numDDDimensions = set_dd_dim(numDomains, ddSettings, &ddGridSetup.ddDimensions);
 953
 954     return ddGridSetup;
 955 }