src/gromacs/tools/tune_pme.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2009-2018, The GROMACS development team.
   5  * Copyright (c) 2019,2020,2021, by the GROMACS development team, led by
   6  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   7  * and including many others, as listed in the AUTHORS file in the
   8  * top-level source directory and at http://www.gromacs.org.
   9  *
  10  * GROMACS is free software; you can redistribute it and/or
  11  * modify it under the terms of the GNU Lesser General Public License
  12  * as published by the Free Software Foundation; either version 2.1
  13  * of the License, or (at your option) any later version.
  14  *
  15  * GROMACS is distributed in the hope that it will be useful,
  16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  18  * Lesser General Public License for more details.
  19  *
  20  * You should have received a copy of the GNU Lesser General Public
  21  * License along with GROMACS; if not, see
  22  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  23  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  24  *
  25  * If you want to redistribute modifications to GROMACS, please
  26  * consider that scientific software is very special. Version
  27  * control is crucial - bugs must be traceable. We will be happy to
  28  * consider code for inclusion in the official distribution, but
  29  * derived work must not be called official GROMACS. Details are found
  30  * in the README & COPYING files - if they are missing, get the
  31  * official version at http://www.gromacs.org.
  32  *
  33  * To help us fund GROMACS development, we humbly ask that you cite
  34  * the research papers on the package. Check out http://www.gromacs.org.
  35  */
  36 #include "gmxpre.h"
  37
  38 #include "tune_pme.h"
  39
  40 #include "config.h"
  41
  42 #include <cmath>
  43 #include <cstdlib>
  44 #include <cstring>
  45 #include <ctime>
  46
  47 #include <algorithm>
  48 #include <numeric>
  49 #include <string>
  50
  51 #ifdef HAVE_SYS_TIME_H
  52 #    include <sys/time.h>
  53 #endif
  54
  55 #include "gromacs/commandline/pargs.h"
  56 #include "gromacs/ewald/pme.h"
  57 #include "gromacs/fft/calcgrid.h"
  58 #include "gromacs/fileio/checkpoint.h"
  59 #include "gromacs/fileio/tpxio.h"
  60 #include "gromacs/math/utilities.h"
  61 #include "gromacs/math/vec.h"
  62 #include "gromacs/mdlib/perf_est.h"
  63 #include "gromacs/mdtypes/commrec.h"
  64 #include "gromacs/mdtypes/inputrec.h"
  65 #include "gromacs/mdtypes/md_enums.h"
  66 #include "gromacs/mdtypes/state.h"
  67 #include "gromacs/pbcutil/pbc.h"
  68 #include "gromacs/taskassignment/usergpuids.h"
  69 #include "gromacs/timing/walltime_accounting.h"
  70 #include "gromacs/topology/topology.h"
  71 #include "gromacs/utility/arraysize.h"
  72 #include "gromacs/utility/baseversion.h"
  73 #include "gromacs/utility/cstringutil.h"
  74 #include "gromacs/utility/fatalerror.h"
  75 #include "gromacs/utility/futil.h"
  76 #include "gromacs/utility/gmxassert.h"
  77 #include "gromacs/utility/path.h"
  78 #include "gromacs/utility/smalloc.h"
  79 #include "gromacs/utility/stringutil.h"
  80
  81 /* Enum for situations that can occur during log file parsing, the
  82  * corresponding string entries can be found in do_the_tests() in
  83  * const char* ParseLog[] */
  84 /* TODO clean up CamelCasing of these enum names */
  85 enum
  86 {
  87     eParselogOK,
  88     eParselogNotFound,
  89     eParselogNoPerfData,
  90     eParselogTerm,
  91     eParselogResetProblem,
  92     eParselogNoDDGrid,
  93     eParselogTPXVersion,
  94     eParselogNotParallel,
  95     eParselogLargePrimeFactor,
  96     eParselogMismatchOfNumberOfPPRanksAndAvailableGPUs,
  97     eParselogGpuProblem,
  98     eParselogFatal,
  99     eParselogNr
 100 };
 101
 102
 103 typedef struct
 104 {
 105     int     nPMEnodes;  /* number of PME-only nodes used in this test */
 106     int     nx, ny, nz; /* DD grid */
 107     int     guessPME;   /* if nPMEnodes == -1, this is the guessed number of PME nodes */
 108     double* Gcycles;    /* This can contain more than one value if doing multiple tests */
 109     double  Gcycles_Av;
 110     float*  ns_per_day;
 111     float   ns_per_day_Av;
 112     float*  PME_f_load;     /* PME mesh/force load average*/
 113     float   PME_f_load_Av;  /* Average average ;) ... */
 114     char*   mdrun_cmd_line; /* Mdrun command line used for this test */
 115 } t_perf;
 116
 117
 118 typedef struct
 119 {
 120     int     nr_inputfiles;  /* The number of tpr and mdp input files */
 121     int64_t orig_sim_steps; /* Number of steps to be done in the real simulation */
 122     int64_t orig_init_step; /* Init step for the real simulation */
 123     real*   rcoulomb;       /* The coulomb radii [0...nr_inputfiles] */
 124     real*   rvdw;           /* The vdW radii */
 125     real*   rlist;          /* Neighbourlist cutoff radius */
 126     int *   nkx, *nky, *nkz;
 127     real *  fsx, *fsy, *fsz; /* Fourierspacing in x,y,z dimension */
 128 } t_inputinfo;
 129
 130
 131 static void sep_line(FILE* fp)
 132 {
 133     fprintf(fp, "\n------------------------------------------------------------\n");
 134 }
 135
 136
 137 /* Wrapper for system calls */
 138 static int gmx_system_call(char* command)
 139 {
 140     return (system(command));
 141 }
 142
 143
 144 /* Check if string starts with substring */
 145 static gmx_bool str_starts(const char* string, const char* substring)
 146 {
 147     return (std::strncmp(string, substring, std::strlen(substring)) == 0);
 148 }
 149
 150
 151 static void cleandata(t_perf* perfdata, int test_nr)
 152 {
 153     perfdata->Gcycles[test_nr]    = 0.0;
 154     perfdata->ns_per_day[test_nr] = 0.0;
 155     perfdata->PME_f_load[test_nr] = 0.0;
 156 }
 157
 158
 159 static void remove_if_exists(const char* fn)
 160 {
 161     if (gmx_fexist(fn))
 162     {
 163         fprintf(stdout, "Deleting %s\n", fn);
 164         remove(fn);
 165     }
 166 }
 167
 168
 169 static void finalize(const char* fn_out)
 170 {
 171     char  buf[STRLEN];
 172     FILE* fp;
 173
 174
 175     fp = fopen(fn_out, "r");
 176     fprintf(stdout, "\n\n");
 177
 178     while (fgets(buf, STRLEN - 1, fp) != nullptr)
 179     {
 180         fprintf(stdout, "%s", buf);
 181     }
 182     fclose(fp);
 183     fprintf(stdout, "\n\n");
 184 }
 185
 186
 187 enum
 188 {
 189     eFoundNothing,
 190     eFoundDDStr,
 191     eFoundAccountingStr,
 192     eFoundCycleStr
 193 };
 194
 195 static int parse_logfile(const char* logfile,
 196                          const char* errfile,
 197                          t_perf*     perfdata,
 198                          int         test_nr,
 199                          int         presteps,
 200                          int64_t     cpt_steps,
 201                          int         nnodes)
 202 {
 203     FILE*      fp;
 204     char       line[STRLEN], dumstring[STRLEN], dumstring2[STRLEN];
 205     const char matchstrdd[]  = "Domain decomposition grid";
 206     const char matchstrcr[]  = "resetting all time and cycle counters";
 207     const char matchstrbal[] = "Average PME mesh/force load:";
 208     const char matchstring[] = "R E A L   C Y C L E   A N D   T I M E   A C C O U N T I N G";
 209     const char errSIG[]      = "signal, stopping at the next";
 210     int        iFound;
 211     float      dum1, dum2, dum3, dum4;
 212     int        ndum;
 213     int        npme;
 214     int64_t    resetsteps     = -1;
 215     gmx_bool   bFoundResetStr = FALSE;
 216     gmx_bool   bResetChecked  = FALSE;
 217
 218
 219     if (!gmx_fexist(logfile))
 220     {
 221         fprintf(stderr, "WARNING: Could not find logfile %s.\n", logfile);
 222         cleandata(perfdata, test_nr);
 223         return eParselogNotFound;
 224     }
 225
 226     fp                            = fopen(logfile, "r");
 227     perfdata->PME_f_load[test_nr] = -1.0;
 228     perfdata->guessPME            = -1;
 229
 230     iFound = eFoundNothing;
 231     if (1 == nnodes)
 232     {
 233         iFound = eFoundDDStr; /* Skip some case statements */
 234     }
 235
 236     while (fgets(line, STRLEN, fp) != nullptr)
 237     {
 238         /* Remove leading spaces */
 239         ltrim(line);
 240
 241         /* Check for TERM and INT signals from user: */
 242         if (std::strstr(line, errSIG) != nullptr)
 243         {
 244             fclose(fp);
 245             cleandata(perfdata, test_nr);
 246             return eParselogTerm;
 247         }
 248
 249         /* Check whether cycle resetting  worked */
 250         if (presteps > 0 && !bFoundResetStr)
 251         {
 252             if (std::strstr(line, matchstrcr) != nullptr)
 253             {
 254                 sprintf(dumstring, "step %s", "%" SCNd64);
 255                 sscanf(line, dumstring, &resetsteps);
 256                 bFoundResetStr = TRUE;
 257                 if (resetsteps == presteps + cpt_steps)
 258                 {
 259                     bResetChecked = TRUE;
 260                 }
 261                 else
 262                 {
 263                     sprintf(dumstring, "%" PRId64, resetsteps);
 264                     sprintf(dumstring2, "%" PRId64, presteps + cpt_steps);
 265                     fprintf(stderr,
 266                             "WARNING: Time step counters were reset at step %s,\n"
 267                             "         though they were supposed to be reset at step %s!\n",
 268                             dumstring,
 269                             dumstring2);
 270                 }
 271             }
 272         }
 273
 274         /* Look for strings that appear in a certain order in the log file: */
 275         switch (iFound)
 276         {
 277             case eFoundNothing:
 278                 /* Look for domain decomp grid and separate PME nodes: */
 279                 if (str_starts(line, matchstrdd))
 280                 {
 281                     sscanf(line,
 282                            "Domain decomposition grid %d x %d x %d, separate PME ranks %d",
 283                            &(perfdata->nx),
 284                            &(perfdata->ny),
 285                            &(perfdata->nz),
 286                            &npme);
 287                     if (perfdata->nPMEnodes == -1)
 288                     {
 289                         perfdata->guessPME = npme;
 290                     }
 291                     else if (perfdata->nPMEnodes != npme)
 292                     {
 293                         gmx_fatal(FARGS,
 294                                   "PME ranks from command line and output file are not identical");
 295                     }
 296                     iFound = eFoundDDStr;
 297                 }
 298                 /* Catch a few errors that might have occurred: */
 299                 else if (str_starts(line, "There is no domain decomposition for"))
 300                 {
 301                     fclose(fp);
 302                     return eParselogNoDDGrid;
 303                 }
 304                 else if (str_starts(line, "The number of ranks you selected"))
 305                 {
 306                     fclose(fp);
 307                     return eParselogLargePrimeFactor;
 308                 }
 309                 else if (str_starts(line, "reading tpx file"))
 310                 {
 311                     fclose(fp);
 312                     return eParselogTPXVersion;
 313                 }
 314                 else if (str_starts(line, "The -dd or -npme option request a parallel simulation"))
 315                 {
 316                     fclose(fp);
 317                     return eParselogNotParallel;
 318                 }
 319                 break;
 320             case eFoundDDStr:
 321                 /* Even after the "Domain decomposition grid" string was found,
 322                  * it could be that mdrun had to quit due to some error. */
 323                 if (str_starts(line, "Incorrect launch configuration: mismatching number of"))
 324                 {
 325                     fclose(fp);
 326                     return eParselogMismatchOfNumberOfPPRanksAndAvailableGPUs;
 327                 }
 328                 else if (str_starts(line, "Some of the requested GPUs do not exist"))
 329                 {
 330                     fclose(fp);
 331                     return eParselogGpuProblem;
 332                 }
 333                 /* Look for PME mesh/force balance (not necessarily present, though) */
 334                 else if (str_starts(line, matchstrbal))
 335                 {
 336                     sscanf(&line[std::strlen(matchstrbal)], "%f", &(perfdata->PME_f_load[test_nr]));
 337                 }
 338                 /* Look for matchstring */
 339                 else if (str_starts(line, matchstring))
 340                 {
 341                     iFound = eFoundAccountingStr;
 342                 }
 343                 break;
 344             case eFoundAccountingStr:
 345                 /* Already found matchstring - look for cycle data */
 346                 if (str_starts(line, "Total  "))
 347                 {
 348                     sscanf(line, "Total %*f %lf", &(perfdata->Gcycles[test_nr]));
 349                     iFound = eFoundCycleStr;
 350                 }
 351                 break;
 352             case eFoundCycleStr:
 353                 /* Already found cycle data - look for remaining performance info and return */
 354                 if (str_starts(line, "Performance:"))
 355                 {
 356                     ndum = sscanf(line, "%s %f %f %f %f", dumstring, &dum1, &dum2, &dum3, &dum4);
 357                     /* (ns/day) is the second last entry, depending on whether GMX_DETAILED_PERF_STATS was set in print_perf(), nrnb.c */
 358                     perfdata->ns_per_day[test_nr] = (ndum == 5) ? dum3 : dum1;
 359                     fclose(fp);
 360                     if (bResetChecked || presteps == 0)
 361                     {
 362                         return eParselogOK;
 363                     }
 364                     else
 365                     {
 366                         return eParselogResetProblem;
 367                     }
 368                 }
 369                 break;
 370         }
 371     } /* while */
 372
 373     /* Close the log file */
 374     fclose(fp);
 375
 376     /* Check why there is no performance data in the log file.
 377      * Did a fatal errors occur? */
 378     if (gmx_fexist(errfile))
 379     {
 380         fp = fopen(errfile, "r");
 381         while (fgets(line, STRLEN, fp) != nullptr)
 382         {
 383             if (str_starts(line, "Fatal error:"))
 384             {
 385                 if (fgets(line, STRLEN, fp) != nullptr)
 386                 {
 387                     fprintf(stderr,
 388                             "\nWARNING: An error occurred during this benchmark:\n"
 389                             "%s\n",
 390                             line);
 391                 }
 392                 fclose(fp);
 393                 cleandata(perfdata, test_nr);
 394                 return eParselogFatal;
 395             }
 396         }
 397         fclose(fp);
 398     }
 399     else
 400     {
 401         fprintf(stderr, "WARNING: Could not find stderr file %s.\n", errfile);
 402     }
 403
 404     /* Giving up ... we could not find out why there is no performance data in
 405      * the log file. */
 406     fprintf(stdout, "No performance data in log file.\n");
 407     cleandata(perfdata, test_nr);
 408
 409     return eParselogNoPerfData;
 410 }
 411
 412
 413 static gmx_bool analyze_data(FILE*        fp,
 414                              const char*  fn,
 415                              t_perf**     perfdata,
 416                              int          nnodes,
 417                              int          ntprs,
 418                              int          ntests,
 419                              int          nrepeats,
 420                              t_inputinfo* info,
 421                              int*         index_tpr, /* OUT: Nr of mdp file with best settings */
 422                              int*         npme_optimal)      /* OUT: Optimal number of PME nodes */
 423 {
 424     int      i, j, k;
 425     int      line = 0, line_win = -1;
 426     int      k_win = -1, i_win = -1, winPME;
 427     double   s = 0.0; /* standard deviation */
 428     t_perf*  pd;
 429     char     strbuf[STRLEN];
 430     char     str_PME_f_load[13];
 431     gmx_bool bCanUseOrigTPR;
 432     gmx_bool bRefinedCoul, bRefinedVdW, bRefinedGrid;
 433
 434
 435     if (nrepeats > 1)
 436     {
 437         sep_line(fp);
 438         fprintf(fp, "Summary of successful runs:\n");
 439         fprintf(fp, "Line tpr PME ranks  Gcycles Av.     Std.dev.       ns/day        PME/f");
 440         if (nnodes > 1)
 441         {
 442             fprintf(fp, "    DD grid");
 443         }
 444         fprintf(fp, "\n");
 445     }
 446
 447
 448     for (k = 0; k < ntprs; k++)
 449     {
 450         for (i = 0; i < ntests; i++)
 451         {
 452             /* Select the right dataset: */
 453             pd = &(perfdata[k][i]);
 454
 455             pd->Gcycles_Av    = 0.0;
 456             pd->PME_f_load_Av = 0.0;
 457             pd->ns_per_day_Av = 0.0;
 458
 459             if (pd->nPMEnodes == -1)
 460             {
 461                 sprintf(strbuf, "(%3d)", pd->guessPME);
 462             }
 463             else
 464             {
 465                 sprintf(strbuf, "     ");
 466             }
 467
 468             /* Get the average run time of a setting */
 469             for (j = 0; j < nrepeats; j++)
 470             {
 471                 pd->Gcycles_Av += pd->Gcycles[j];
 472                 pd->PME_f_load_Av += pd->PME_f_load[j];
 473             }
 474             pd->Gcycles_Av /= nrepeats;
 475             pd->PME_f_load_Av /= nrepeats;
 476
 477             for (j = 0; j < nrepeats; j++)
 478             {
 479                 if (pd->ns_per_day[j] > 0.0)
 480                 {
 481                     pd->ns_per_day_Av += pd->ns_per_day[j];
 482                 }
 483                 else
 484                 {
 485                     /* Somehow the performance number was not aquired for this run,
 486                      * therefor set the average to some negative value: */
 487                     pd->ns_per_day_Av = -1.0F * nrepeats;
 488                     break;
 489                 }
 490             }
 491             pd->ns_per_day_Av /= nrepeats;
 492
 493             /* Nicer output: */
 494             if (pd->PME_f_load_Av > 0.0)
 495             {
 496                 sprintf(str_PME_f_load, "%12.3f", pd->PME_f_load_Av);
 497             }
 498             else
 499             {
 500                 sprintf(str_PME_f_load, "%s", "         -  ");
 501             }
 502
 503
 504             /* We assume we had a successful run if both averages are positive */
 505             if (pd->Gcycles_Av > 0.0 && pd->ns_per_day_Av > 0.0)
 506             {
 507                 /* Output statistics if repeats were done */
 508                 if (nrepeats > 1)
 509                 {
 510                     /* Calculate the standard deviation */
 511                     s = 0.0;
 512                     for (j = 0; j < nrepeats; j++)
 513                     {
 514                         s += gmx::square(pd->Gcycles[j] - pd->Gcycles_Av);
 515                     }
 516                     s /= (nrepeats - 1);
 517                     s = std::sqrt(s);
 518
 519                     fprintf(fp,
 520                             "%4d %3d %4d%s %12.3f %12.3f %12.3f %s",
 521                             line,
 522                             k,
 523                             pd->nPMEnodes,
 524                             strbuf,
 525                             pd->Gcycles_Av,
 526                             s,
 527                             pd->ns_per_day_Av,
 528                             str_PME_f_load);
 529                     if (nnodes > 1)
 530                     {
 531                         fprintf(fp, "  %3d %3d %3d", pd->nx, pd->ny, pd->nz);
 532                     }
 533                     fprintf(fp, "\n");
 534                 }
 535                 /* Store the index of the best run found so far in 'winner': */
 536                 if ((k_win == -1) || (pd->Gcycles_Av < perfdata[k_win][i_win].Gcycles_Av))
 537                 {
 538                     k_win    = k;
 539                     i_win    = i;
 540                     line_win = line;
 541                 }
 542                 line++;
 543             }
 544         }
 545     }
 546
 547     if (k_win == -1)
 548     {
 549         gmx_fatal(FARGS, "None of the runs was successful! Check %s for problems.", fn);
 550     }
 551
 552     sep_line(fp);
 553
 554     winPME = perfdata[k_win][i_win].nPMEnodes;
 555
 556     if (1 == ntests)
 557     {
 558         /* We stuck to a fixed number of PME-only nodes */
 559         sprintf(strbuf, "settings No. %d", k_win);
 560     }
 561     else
 562     {
 563         /* We have optimized the number of PME-only nodes */
 564         if (winPME == -1)
 565         {
 566             sprintf(strbuf, "%s", "the automatic number of PME ranks");
 567         }
 568         else
 569         {
 570             sprintf(strbuf, "%d PME ranks", winPME);
 571         }
 572     }
 573     fprintf(fp, "Best performance was achieved with %s", strbuf);
 574     if ((nrepeats > 1) && (ntests > 1))
 575     {
 576         fprintf(fp, " (see line %d)", line_win);
 577     }
 578     fprintf(fp, "\n");
 579
 580     /* Only mention settings if they were modified: */
 581     bRefinedCoul = !gmx_within_tol(info->rcoulomb[k_win], info->rcoulomb[0], GMX_REAL_EPS);
 582     bRefinedVdW  = !gmx_within_tol(info->rvdw[k_win], info->rvdw[0], GMX_REAL_EPS);
 583     bRefinedGrid = !(info->nkx[k_win] == info->nkx[0] && info->nky[k_win] == info->nky[0]
 584                      && info->nkz[k_win] == info->nkz[0]);
 585
 586     if (bRefinedCoul || bRefinedVdW || bRefinedGrid)
 587     {
 588         fprintf(fp, "Optimized PME settings:\n");
 589         bCanUseOrigTPR = FALSE;
 590     }
 591     else
 592     {
 593         bCanUseOrigTPR = TRUE;
 594     }
 595
 596     if (bRefinedCoul)
 597     {
 598         fprintf(fp, "   New Coulomb radius: %f nm (was %f nm)\n", info->rcoulomb[k_win], info->rcoulomb[0]);
 599     }
 600
 601     if (bRefinedVdW)
 602     {
 603         fprintf(fp, "   New Van der Waals radius: %f nm (was %f nm)\n", info->rvdw[k_win], info->rvdw[0]);
 604     }
 605
 606     if (bRefinedGrid)
 607     {
 608         fprintf(fp,
 609                 "   New Fourier grid xyz: %d %d %d (was %d %d %d)\n",
 610                 info->nkx[k_win],
 611                 info->nky[k_win],
 612                 info->nkz[k_win],
 613                 info->nkx[0],
 614                 info->nky[0],
 615                 info->nkz[0]);
 616     }
 617
 618     if (bCanUseOrigTPR && ntprs > 1)
 619     {
 620         fprintf(fp, "and original PME settings.\n");
 621     }
 622
 623     fflush(fp);
 624
 625     /* Return the index of the mdp file that showed the highest performance
 626      * and the optimal number of PME nodes */
 627     *index_tpr    = k_win;
 628     *npme_optimal = winPME;
 629
 630     return bCanUseOrigTPR;
 631 }
 632
 633
 634 /* Get the commands we need to set up the runs from environment variables */
 635 static void get_program_paths(gmx_bool bThreads, char* cmd_mpirun[], char* cmd_mdrun[])
 636 {
 637     char*      cp;
 638     const char def_mpirun[] = "mpirun";
 639
 640     const char empty_mpirun[] = "";
 641
 642     /* Get the commands we need to set up the runs from environment variables */
 643     if (!bThreads)
 644     {
 645         if ((cp = getenv("MPIRUN")) != nullptr)
 646         {
 647             *cmd_mpirun = gmx_strdup(cp);
 648         }
 649         else
 650         {
 651             *cmd_mpirun = gmx_strdup(def_mpirun);
 652         }
 653     }
 654     else
 655     {
 656         *cmd_mpirun = gmx_strdup(empty_mpirun);
 657     }
 658
 659     if (*cmd_mdrun == nullptr)
 660     {
 661         /* The use of MDRUN is deprecated, but made available in 5.1
 662            for backward compatibility. It may be removed in a future
 663            version. */
 664         if ((cp = getenv("MDRUN")) != nullptr)
 665         {
 666             *cmd_mdrun = gmx_strdup(cp);
 667         }
 668         else
 669         {
 670             gmx_fatal(FARGS, "The way to call mdrun must be set in the -mdrun command-line flag.");
 671         }
 672     }
 673 }
 674
 675 /* Check that the commands will run mdrun (perhaps via mpirun) by
 676  * running a very quick test simulation. Requires MPI environment or
 677  * GPU support to be available if applicable. */
 678 /* TODO implement feature to parse the log file to get the list of
 679    compatible GPUs from mdrun, if the user of gmx tune-pme has not
 680    given one. */
 681 static void check_mdrun_works(gmx_bool    bThreads,
 682                               const char* cmd_mpirun,
 683                               const char* cmd_np,
 684                               const char* cmd_mdrun,
 685                               gmx_bool    bNeedGpuSupport)
 686 {
 687     char*      command = nullptr;
 688     char*      cp;
 689     char       line[STRLEN];
 690     FILE*      fp;
 691     const char filename[] = "benchtest.log";
 692
 693     /* This string should always be identical to the one in copyrite.c,
 694      * gmx_print_version_info() in the GMX_MPI section */
 695     const char match_mpi[]     = "MPI library:        MPI";
 696     const char match_mdrun[]   = "Executable: ";
 697     const char match_nogpu[]   = "GPU support:        disabled";
 698     gmx_bool   bMdrun          = FALSE;
 699     gmx_bool   bMPI            = FALSE;
 700     gmx_bool   bHaveGpuSupport = TRUE;
 701
 702     /* Run a small test to see whether mpirun + mdrun work  */
 703     fprintf(stdout, "Making sure that mdrun can be executed. ");
 704     if (bThreads)
 705     {
 706         snew(command, std::strlen(cmd_mdrun) + std::strlen(cmd_np) + std::strlen(filename) + 50);
 707         sprintf(command, "%s%s -version -maxh 0.001 1> %s 2>&1", cmd_mdrun, cmd_np, filename);
 708     }
 709     else
 710     {
 711         snew(command,
 712              std::strlen(cmd_mpirun) + std::strlen(cmd_np) + std::strlen(cmd_mdrun)
 713                      + std::strlen(filename) + 50);
 714         sprintf(command, "%s%s%s -version -maxh 0.001 1> %s 2>&1", cmd_mpirun, cmd_np, cmd_mdrun, filename);
 715     }
 716     fprintf(stdout, "Trying '%s' ... ", command);
 717     make_backup(filename);
 718     gmx_system_call(command);
 719
 720     /* Check if we find the characteristic string in the output: */
 721     if (!gmx_fexist(filename))
 722     {
 723         gmx_fatal(FARGS, "Output from test run could not be found.");
 724     }
 725
 726     fp = fopen(filename, "r");
 727     /* We need to scan the whole output file, since sometimes the queuing system
 728      * also writes stuff to stdout/err */
 729     while (!feof(fp))
 730     {
 731         cp = fgets(line, STRLEN, fp);
 732         if (cp != nullptr)
 733         {
 734             if (str_starts(line, match_mdrun))
 735             {
 736                 bMdrun = TRUE;
 737             }
 738             if (str_starts(line, match_mpi))
 739             {
 740                 bMPI = TRUE;
 741             }
 742             if (str_starts(line, match_nogpu))
 743             {
 744                 bHaveGpuSupport = FALSE;
 745             }
 746         }
 747     }
 748     fclose(fp);
 749
 750     if (bThreads)
 751     {
 752         if (bMPI)
 753         {
 754             gmx_fatal(FARGS,
 755                       "Need a threaded version of mdrun. This one\n"
 756                       "(%s)\n"
 757                       "seems to have been compiled with MPI instead.",
 758                       cmd_mdrun);
 759         }
 760     }
 761     else
 762     {
 763         if (bMdrun && !bMPI)
 764         {
 765             gmx_fatal(FARGS,
 766                       "Need an MPI-enabled version of mdrun. This one\n"
 767                       "(%s)\n"
 768                       "seems to have been compiled without MPI support.",
 769                       cmd_mdrun);
 770         }
 771     }
 772
 773     if (!bMdrun)
 774     {
 775         gmx_fatal(FARGS, "Cannot execute mdrun. Please check %s for problems!", filename);
 776     }
 777
 778     if (bNeedGpuSupport && !bHaveGpuSupport)
 779     {
 780         gmx_fatal(FARGS, "The mdrun executable did not have the expected GPU support.");
 781     }
 782
 783     fprintf(stdout, "passed.\n");
 784
 785     /* Clean up ... */
 786     remove(filename);
 787     sfree(command);
 788 }
 789
 790 /* Handles the no-GPU case by emitting an empty string. */
 791 static std::string make_gpu_id_command_line(const char* eligible_gpu_ids)
 792 {
 793     /* If the user has given no eligible GPU IDs, or we're trying the
 794      * default behaviour, then there is nothing for tune_pme to give
 795      * to mdrun -gpu_id */
 796     if (eligible_gpu_ids != nullptr)
 797     {
 798         return gmx::formatString("-gpu_id %s", eligible_gpu_ids);
 799     }
 800
 801
 802     return std::string();
 803 }
 804
 805 static void launch_simulation(gmx_bool    bLaunch,          /* Should the simulation be launched? */
 806                               FILE*       fp,               /* General log file */
 807                               gmx_bool    bThreads,         /* whether to use threads */
 808                               char*       cmd_mpirun,       /* Command for mpirun */
 809                               char*       cmd_np,           /* Switch for -np or -ntmpi or empty */
 810                               char*       cmd_mdrun,        /* Command for mdrun */
 811                               char*       args_for_mdrun,   /* Arguments for mdrun */
 812                               const char* simulation_tpr,   /* This tpr will be simulated */
 813                               int         nPMEnodes,        /* Number of PME ranks to use */
 814                               const char* eligible_gpu_ids) /* Available GPU IDs for
 815                                                              * constructing mdrun command lines */
 816 {
 817     char* command;
 818
 819
 820     /* Make enough space for the system call command,
 821      * (200 extra chars for -npme ... etc. options should suffice): */
 822     snew(command,
 823          std::strlen(cmd_mpirun) + std::strlen(cmd_mdrun) + std::strlen(cmd_np)
 824                  + std::strlen(args_for_mdrun) + std::strlen(simulation_tpr) + 200);
 825
 826     auto cmd_gpu_ids = make_gpu_id_command_line(eligible_gpu_ids);
 827
 828     /* Note that the -passall options requires args_for_mdrun to be at the end
 829      * of the command line string */
 830     if (bThreads)
 831     {
 832         sprintf(command,
 833                 "%s%s-npme %d -s %s %s %s",
 834                 cmd_mdrun,
 835                 cmd_np,
 836                 nPMEnodes,
 837                 simulation_tpr,
 838                 args_for_mdrun,
 839                 cmd_gpu_ids.c_str());
 840     }
 841     else
 842     {
 843         sprintf(command,
 844                 "%s%s%s -npme %d -s %s %s %s",
 845                 cmd_mpirun,
 846                 cmd_np,
 847                 cmd_mdrun,
 848                 nPMEnodes,
 849                 simulation_tpr,
 850                 args_for_mdrun,
 851                 cmd_gpu_ids.c_str());
 852     }
 853
 854     fprintf(fp, "%s this command line to launch the simulation:\n\n%s", bLaunch ? "Using" : "Please use", command);
 855     sep_line(fp);
 856     fflush(fp);
 857
 858     /* Now the real thing! */
 859     if (bLaunch)
 860     {
 861         fprintf(stdout, "\nLaunching simulation with best parameters now.\nExecuting '%s'", command);
 862         sep_line(stdout);
 863         fflush(stdout);
 864         gmx_system_call(command);
 865     }
 866 }
 867
 868
 869 static void modify_PMEsettings(int64_t     simsteps,    /* Set this value as number of time steps */
 870                                int64_t     init_step,   /* Set this value as init_step */
 871                                const char* fn_best_tpr, /* tpr file with the best performance */
 872                                const char* fn_sim_tpr)  /* name of tpr file to be launched */
 873 {
 874     t_state    state;
 875     gmx_mtop_t mtop;
 876     char       buf[200];
 877
 878     t_inputrec  irInstance;
 879     t_inputrec* ir = &irInstance;
 880     read_tpx_state(fn_best_tpr, ir, &state, &mtop);
 881
 882     /* Reset nsteps and init_step to the value of the input .tpr file */
 883     ir->nsteps    = simsteps;
 884     ir->init_step = init_step;
 885
 886     /* Write the tpr file which will be launched */
 887     sprintf(buf, "Writing optimized simulation file %s with nsteps=%s.\n", fn_sim_tpr, "%" PRId64);
 888     fprintf(stdout, buf, ir->nsteps);
 889     fflush(stdout);
 890     write_tpx_state(fn_sim_tpr, ir, &state, mtop);
 891 }
 892
 893 static gmx_bool can_scale_rvdw(VanDerWaalsType vdwtype)
 894 {
 895     return (VanDerWaalsType::Cut == vdwtype || VanDerWaalsType::Pme == vdwtype);
 896 }
 897
 898 #define EPME_SWITCHED(e) \
 899     ((e) == CoulombInteractionType::PmeSwitch || (e) == CoulombInteractionType::PmeUserSwitch)
 900
 901 /* Make additional TPR files with more computational load for the
 902  * direct space processors: */
 903 static void make_benchmark_tprs(const char* fn_sim_tpr,  /* READ : User-provided tpr file */
 904                                 char*   fn_bench_tprs[], /* WRITE: Names of benchmark tpr files  */
 905                                 int64_t benchsteps, /* Number of time steps for benchmark runs */
 906                                 int64_t statesteps, /* Step counter in checkpoint file */
 907                                 real    rmin,    /* Minimal Coulomb radius                        */
 908                                 real    rmax,    /* Maximal Coulomb radius                        */
 909                                 bool bScaleRvdw, /* Scale rvdw along with rcoulomb                */
 910                                 const int* ntprs,  /* No. of TPRs to write, each with a different
 911                                                       rcoulomb and fourierspacing  */
 912                                 t_inputinfo* info, /* Contains information about mdp file options */
 913                                 FILE*        fp) /* Write the output here                         */
 914 {
 915     int        i, j, d;
 916     t_state    state;
 917     gmx_mtop_t mtop;
 918     real       nlist_buffer; /* Thickness of the buffer regions for PME-switch potentials */
 919     char       buf[200];
 920     rvec       box_size;
 921     gmx_bool   bNote = FALSE;
 922     real       add;            /* Add this to rcoul for the next test    */
 923     real       fac = 1.0;      /* Scaling factor for Coulomb radius      */
 924     real       fourierspacing; /* Basic fourierspacing from tpr          */
 925
 926
 927     sprintf(buf,
 928             "Making benchmark tpr file%s with %s time step%s",
 929             *ntprs > 1 ? "s" : "",
 930             "%" PRId64,
 931             benchsteps > 1 ? "s" : "");
 932     fprintf(stdout, buf, benchsteps);
 933     if (statesteps > 0)
 934     {
 935         sprintf(buf, " (adding %s steps from checkpoint file)", "%" PRId64);
 936         fprintf(stdout, buf, statesteps);
 937         benchsteps += statesteps;
 938     }
 939     fprintf(stdout, ".\n");
 940
 941     t_inputrec  irInstance;
 942     t_inputrec* ir = &irInstance;
 943     read_tpx_state(fn_sim_tpr, ir, &state, &mtop);
 944
 945     /* Check if some kind of PME was chosen */
 946     if (EEL_PME(ir->coulombtype) == FALSE)
 947     {
 948         gmx_fatal(FARGS,
 949                   "Can only do optimizations for simulations with %s electrostatics.",
 950                   enumValueToString(CoulombInteractionType::Pme));
 951     }
 952
 953     /* Check if rcoulomb == rlist, which is necessary for plain PME. */
 954     if ((ir->cutoff_scheme != CutoffScheme::Verlet)
 955         && (CoulombInteractionType::Pme == ir->coulombtype) && !(ir->rcoulomb == ir->rlist))
 956     {
 957         gmx_fatal(FARGS,
 958                   "%s requires rcoulomb (%f) to be equal to rlist (%f).",
 959                   enumValueToString(CoulombInteractionType::Pme),
 960                   ir->rcoulomb,
 961                   ir->rlist);
 962     }
 963     /* For other PME types, rcoulomb is allowed to be smaller than rlist */
 964     else if (ir->rcoulomb > ir->rlist)
 965     {
 966         gmx_fatal(FARGS,
 967                   "%s requires rcoulomb (%f) to be equal to or smaller than rlist (%f)",
 968                   enumValueToString(ir->coulombtype),
 969                   ir->rcoulomb,
 970                   ir->rlist);
 971     }
 972
 973     if (bScaleRvdw && ir->rvdw != ir->rcoulomb)
 974     {
 975         fprintf(stdout, "NOTE: input rvdw != rcoulomb, will not scale rvdw\n");
 976         bScaleRvdw = FALSE;
 977     }
 978
 979     /* Reduce the number of steps for the benchmarks */
 980     info->orig_sim_steps = ir->nsteps;
 981     ir->nsteps           = benchsteps;
 982     /* We must not use init_step from the input tpr file for the benchmarks */
 983     info->orig_init_step = ir->init_step;
 984     ir->init_step        = 0;
 985
 986     /* For PME-switch potentials, keep the radial distance of the buffer region */
 987     nlist_buffer = ir->rlist - ir->rcoulomb;
 988
 989     /* Determine length of triclinic box vectors */
 990     for (d = 0; d < DIM; d++)
 991     {
 992         box_size[d] = 0;
 993         for (i = 0; i < DIM; i++)
 994         {
 995             box_size[d] += state.box[d][i] * state.box[d][i];
 996         }
 997         box_size[d] = std::sqrt(box_size[d]);
 998     }
 999
1000     if (ir->fourier_spacing > 0)
1001     {
1002         info->fsx[0] = ir->fourier_spacing;
1003         info->fsy[0] = ir->fourier_spacing;
1004         info->fsz[0] = ir->fourier_spacing;
1005     }
1006     else
1007     {
1008         /* Reconstruct fourierspacing per dimension from the number of grid points and box size */
1009         info->fsx[0] = box_size[XX] / ir->nkx;
1010         info->fsy[0] = box_size[YY] / ir->nky;
1011         info->fsz[0] = box_size[ZZ] / ir->nkz;
1012     }
1013
1014     /* If no value for the fourierspacing was provided on the command line, we
1015      * use the reconstruction from the tpr file */
1016     if (ir->fourier_spacing > 0)
1017     {
1018         /* Use the spacing from the tpr */
1019         fourierspacing = ir->fourier_spacing;
1020     }
1021     else
1022     {
1023         /* Use the maximum observed spacing */
1024         fourierspacing = std::max(std::max(info->fsx[0], info->fsy[0]), info->fsz[0]);
1025     }
1026
1027     fprintf(stdout, "Calculating PME grid points on the basis of a fourierspacing of %f nm\n", fourierspacing);
1028
1029     /* For performance comparisons the number of particles is useful to have */
1030     fprintf(fp, "   Number of particles  : %d\n", mtop.natoms);
1031
1032     /* Print information about settings of which some are potentially modified: */
1033     fprintf(fp, "   Coulomb type         : %s\n", enumValueToString(ir->coulombtype));
1034     fprintf(fp,
1035             "   Grid spacing x y z   : %f %f %f\n",
1036             box_size[XX] / ir->nkx,
1037             box_size[YY] / ir->nky,
1038             box_size[ZZ] / ir->nkz);
1039     fprintf(fp, "   Van der Waals type   : %s\n", enumValueToString(ir->vdwtype));
1040     if (ir_vdw_switched(ir))
1041     {
1042         fprintf(fp, "   rvdw_switch          : %f nm\n", ir->rvdw_switch);
1043     }
1044     if (EPME_SWITCHED(ir->coulombtype))
1045     {
1046         fprintf(fp, "   rlist                : %f nm\n", ir->rlist);
1047     }
1048
1049     /* Print a descriptive line about the tpr settings tested */
1050     fprintf(fp, "\nWill try these real/reciprocal workload settings:\n");
1051     fprintf(fp, " No.   scaling  rcoulomb");
1052     fprintf(fp, "  nkx  nky  nkz");
1053     fprintf(fp, "   spacing");
1054     if (can_scale_rvdw(ir->vdwtype))
1055     {
1056         fprintf(fp, "      rvdw");
1057     }
1058     if (EPME_SWITCHED(ir->coulombtype))
1059     {
1060         fprintf(fp, "     rlist");
1061     }
1062     fprintf(fp, "  tpr file\n");
1063
1064     /* Loop to create the requested number of tpr input files */
1065     for (j = 0; j < *ntprs; j++)
1066     {
1067         /* The first .tpr is the provided one, just need to modify nsteps,
1068          * so skip the following block */
1069         if (j != 0)
1070         {
1071             /* Determine which Coulomb radii rc to use in the benchmarks */
1072             add = (rmax - rmin) / (*ntprs - 1);
1073             if (gmx_within_tol(rmin, info->rcoulomb[0], GMX_REAL_EPS))
1074             {
1075                 ir->rcoulomb = rmin + j * add;
1076             }
1077             else if (gmx_within_tol(rmax, info->rcoulomb[0], GMX_REAL_EPS))
1078             {
1079                 ir->rcoulomb = rmin + (j - 1) * add;
1080             }
1081             else
1082             {
1083                 /* rmin != rcoul != rmax, ergo test between rmin and rmax */
1084                 add          = (rmax - rmin) / (*ntprs - 2);
1085                 ir->rcoulomb = rmin + (j - 1) * add;
1086             }
1087
1088             /* Determine the scaling factor fac */
1089             fac = ir->rcoulomb / info->rcoulomb[0];
1090
1091             /* Scale the Fourier grid spacing */
1092             ir->nkx = ir->nky = ir->nkz = 0;
1093             calcFftGrid(nullptr,
1094                         state.box,
1095                         fourierspacing * fac,
1096                         minimalPmeGridSize(ir->pme_order),
1097                         &ir->nkx,
1098                         &ir->nky,
1099                         &ir->nkz);
1100
1101             /* Adjust other radii since various conditions need to be fulfilled */
1102             if (CoulombInteractionType::Pme == ir->coulombtype)
1103             {
1104                 /* plain PME, rcoulomb must be equal to rlist TODO only in the group scheme? */
1105                 ir->rlist = ir->rcoulomb;
1106             }
1107             else
1108             {
1109                 /* rlist must be >= rcoulomb, we keep the size of the buffer region */
1110                 ir->rlist = ir->rcoulomb + nlist_buffer;
1111             }
1112
1113             if (bScaleRvdw && can_scale_rvdw(ir->vdwtype))
1114             {
1115                 if (CutoffScheme::Verlet == ir->cutoff_scheme || VanDerWaalsType::Pme == ir->vdwtype)
1116                 {
1117                     /* With either the Verlet cutoff-scheme or LJ-PME,
1118                        the van der Waals radius must always equal the
1119                        Coulomb radius */
1120                     ir->rvdw = ir->rcoulomb;
1121                 }
1122                 else
1123                 {
1124                     /* For vdw cutoff, rvdw >= rlist */
1125                     ir->rvdw = std::max(info->rvdw[0], ir->rlist);
1126                 }
1127             }
1128         } /* end of "if (j != 0)" */
1129
1130         /* for j==0: Save the original settings
1131          * for j >0: Save modified radii and Fourier grids */
1132         info->rcoulomb[j] = ir->rcoulomb;
1133         info->rvdw[j]     = ir->rvdw;
1134         info->nkx[j]      = ir->nkx;
1135         info->nky[j]      = ir->nky;
1136         info->nkz[j]      = ir->nkz;
1137         info->rlist[j]    = ir->rlist;
1138         info->fsx[j]      = fac * fourierspacing;
1139         info->fsy[j]      = fac * fourierspacing;
1140         info->fsz[j]      = fac * fourierspacing;
1141
1142         /* Write the benchmark tpr file */
1143         fn_bench_tprs[j] = gmx_strdup(
1144                 gmx::Path::concatenateBeforeExtension(fn_sim_tpr, gmx::formatString("_bench%.2d", j))
1145                         .c_str());
1146
1147         fprintf(stdout, "Writing benchmark tpr %s with nsteps=", fn_bench_tprs[j]);
1148         fprintf(stdout, "%" PRId64, ir->nsteps);
1149         if (j > 0)
1150         {
1151             fprintf(stdout, ", scaling factor %f\n", fac);
1152         }
1153         else
1154         {
1155             fprintf(stdout, ", unmodified settings\n");
1156         }
1157
1158         write_tpx_state(fn_bench_tprs[j], ir, &state, mtop);
1159
1160         /* Write information about modified tpr settings to log file */
1161         fprintf(fp, "%4d%10f%10f", j, fac, ir->rcoulomb);
1162         fprintf(fp, "%5d%5d%5d", ir->nkx, ir->nky, ir->nkz);
1163         fprintf(fp, " %9f ", info->fsx[j]);
1164         if (can_scale_rvdw(ir->vdwtype))
1165         {
1166             fprintf(fp, "%10f", ir->rvdw);
1167         }
1168         if (EPME_SWITCHED(ir->coulombtype))
1169         {
1170             fprintf(fp, "%10f", ir->rlist);
1171         }
1172         fprintf(fp, "  %-14s\n", fn_bench_tprs[j]);
1173
1174         /* Make it clear to the user that some additional settings were modified */
1175         if (!gmx_within_tol(ir->rvdw, info->rvdw[0], GMX_REAL_EPS)
1176             || !gmx_within_tol(ir->rlist, info->rlist[0], GMX_REAL_EPS))
1177         {
1178             bNote = TRUE;
1179         }
1180     }
1181     if (bNote)
1182     {
1183         fprintf(fp,
1184                 "\nNote that in addition to the Coulomb radius and the Fourier grid\n"
1185                 "other input settings were also changed (see table above).\n"
1186                 "Please check if the modified settings are appropriate.\n");
1187     }
1188     fflush(stdout);
1189     fflush(fp);
1190 }
1191
1192
1193 /* Rename the files we want to keep to some meaningful filename and
1194  * delete the rest */
1195 static void cleanup(const t_filenm* fnm, int nfile, int k, int nnodes, int nPMEnodes, int nr, gmx_bool bKeepStderr)
1196 {
1197     char        numstring[STRLEN];
1198     const char* fn = nullptr;
1199     int         i;
1200     const char* opt;
1201
1202
1203     fprintf(stdout, "Cleaning up, deleting benchmark temp files ...\n");
1204
1205     for (i = 0; i < nfile; i++)
1206     {
1207         opt = const_cast<char*>(fnm[i].opt);
1208         if (std::strcmp(opt, "-p") == 0)
1209         {
1210             /* do nothing; keep this file */
1211         }
1212         else if (std::strcmp(opt, "-bg") == 0)
1213         {
1214             /* Give the log file a nice name so one can later see which parameters were used */
1215             numstring[0] = '\0';
1216             if (nr > 0)
1217             {
1218                 sprintf(numstring, "_%d", nr);
1219             }
1220             std::string newfilename = gmx::formatString(
1221                     "%s_no%d_np%d_npme%d%s", opt2fn("-bg", nfile, fnm), k, nnodes, nPMEnodes, numstring);
1222             if (gmx_fexist(opt2fn("-bg", nfile, fnm)))
1223             {
1224                 fprintf(stdout, "renaming log file to %s\n", newfilename.c_str());
1225                 make_backup(newfilename);
1226                 rename(opt2fn("-bg", nfile, fnm), newfilename.c_str());
1227             }
1228         }
1229         else if (std::strcmp(opt, "-err") == 0)
1230         {
1231             /* This file contains the output of stderr. We want to keep it in
1232              * cases where there have been problems. */
1233             fn           = opt2fn(opt, nfile, fnm);
1234             numstring[0] = '\0';
1235             if (nr > 0)
1236             {
1237                 sprintf(numstring, "_%d", nr);
1238             }
1239             std::string newfilename =
1240                     gmx::formatString("%s_no%d_np%d_npme%d%s", fn, k, nnodes, nPMEnodes, numstring);
1241             if (gmx_fexist(fn))
1242             {
1243                 if (bKeepStderr)
1244                 {
1245                     fprintf(stdout, "Saving stderr output in %s\n", newfilename.c_str());
1246                     make_backup(newfilename);
1247                     rename(fn, newfilename.c_str());
1248                 }
1249                 else
1250                 {
1251                     fprintf(stdout, "Deleting %s\n", fn);
1252                     remove(fn);
1253                 }
1254             }
1255         }
1256         /* Delete the files which are created for each benchmark run: (options -b*) */
1257         else if ((0 == std::strncmp(opt, "-b", 2)) && (opt2bSet(opt, nfile, fnm) || !is_optional(&fnm[i])))
1258         {
1259             remove_if_exists(opt2fn(opt, nfile, fnm));
1260         }
1261     }
1262 }
1263
1264
1265 enum
1266 {
1267     eNpmeAuto,
1268     eNpmeAll,
1269     eNpmeReduced,
1270     eNpmeSubset,
1271     eNpmeNr
1272 };
1273
1274 /* Create a list of numbers of PME nodes to test */
1275 static void make_npme_list(const char* npmevalues_opt, /* Make a complete list with all
1276                                                         * possibilities or a short list that keeps
1277                                                         * only reasonable numbers of PME nodes */
1278                            int* nentries,    /* Number of entries we put in the nPMEnodes list   */
1279                            int* nPMEnodes[], /* Each entry contains the value for -npme          */
1280                            int  nnodes,      /* Total number of nodes to do the tests on         */
1281                            int  minPMEnodes, /* Minimum number of PME nodes                      */
1282                            int  maxPMEnodes)  /* Maximum number of PME nodes                      */
1283 {
1284     int i, npme, npp;
1285     int min_factor = 1; /* We request that npp and npme have this minimal
1286                          * largest common factor (depends on npp)           */
1287     int nlistmax;       /* Max. list size                                   */
1288     int nlist;          /* Actual number of entries in list                 */
1289     int eNPME = 0;
1290
1291
1292     /* Do we need to check all possible values for -npme or is a reduced list enough? */
1293     if (!std::strcmp(npmevalues_opt, "all"))
1294     {
1295         eNPME = eNpmeAll;
1296     }
1297     else if (!std::strcmp(npmevalues_opt, "subset"))
1298     {
1299         eNPME = eNpmeSubset;
1300     }
1301     else /* "auto" or "range" */
1302     {
1303         if (nnodes <= 64)
1304         {
1305             eNPME = eNpmeAll;
1306         }
1307         else if (nnodes < 128)
1308         {
1309             eNPME = eNpmeReduced;
1310         }
1311         else
1312         {
1313             eNPME = eNpmeSubset;
1314         }
1315     }
1316
1317     /* Calculate how many entries we could possibly have (in case of -npme all) */
1318     if (nnodes > 2)
1319     {
1320         nlistmax = maxPMEnodes - minPMEnodes + 3;
1321         if (0 == minPMEnodes)
1322         {
1323             nlistmax--;
1324         }
1325     }
1326     else
1327     {
1328         nlistmax = 1;
1329     }
1330
1331     /* Now make the actual list which is at most of size nlist */
1332     snew(*nPMEnodes, nlistmax);
1333     nlist = 0; /* start counting again, now the real entries in the list */
1334     for (i = 0; i < nlistmax - 2; i++)
1335     {
1336         npme = maxPMEnodes - i;
1337         npp  = nnodes - npme;
1338         switch (eNPME)
1339         {
1340             case eNpmeAll: min_factor = 1; break;
1341             case eNpmeReduced: min_factor = 2; break;
1342             case eNpmeSubset:
1343                 /* For 2d PME we want a common largest factor of at least the cube
1344                  * root of the number of PP nodes */
1345                 min_factor = static_cast<int>(std::cbrt(npp));
1346                 break;
1347             default: gmx_fatal(FARGS, "Unknown option for eNPME in make_npme_list");
1348         }
1349         if (std::gcd(npp, npme) >= min_factor)
1350         {
1351             (*nPMEnodes)[nlist] = npme;
1352             nlist++;
1353         }
1354     }
1355     /* We always test 0 PME nodes and the automatic number */
1356     *nentries               = nlist + 2;
1357     (*nPMEnodes)[nlist]     = 0;
1358     (*nPMEnodes)[nlist + 1] = -1;
1359
1360     fprintf(stderr, "Will try the following %d different values for -npme:\n", *nentries);
1361     for (i = 0; i < *nentries - 1; i++)
1362     {
1363         fprintf(stderr, "%d, ", (*nPMEnodes)[i]);
1364     }
1365     fprintf(stderr, "and %d (auto).\n", (*nPMEnodes)[*nentries - 1]);
1366 }
1367
1368
1369 /* Allocate memory to store the performance data */
1370 static void init_perfdata(t_perf* perfdata[], int ntprs, int datasets, int repeats)
1371 {
1372     int i, j, k;
1373
1374
1375     for (k = 0; k < ntprs; k++)
1376     {
1377         snew(perfdata[k], datasets);
1378         for (i = 0; i < datasets; i++)
1379         {
1380             for (j = 0; j < repeats; j++)
1381             {
1382                 snew(perfdata[k][i].Gcycles, repeats);
1383                 snew(perfdata[k][i].ns_per_day, repeats);
1384                 snew(perfdata[k][i].PME_f_load, repeats);
1385             }
1386         }
1387     }
1388 }
1389
1390
1391 /* Check for errors on mdrun -h */
1392 static void make_sure_it_runs(char* mdrun_cmd_line, int length, FILE* fp, const t_filenm* fnm, int nfile)
1393 {
1394     char *command, *msg;
1395     int   ret;
1396
1397     snew(command, length + 15);
1398     snew(msg, length + 500);
1399
1400     fprintf(stdout, "Making sure the benchmarks can be executed by running just 1 step...\n");
1401     sprintf(command, "%s -nsteps 1 -quiet", mdrun_cmd_line);
1402     fprintf(stdout, "Executing '%s' ...\n", command);
1403     ret = gmx_system_call(command);
1404
1405     if (0 != ret)
1406     {
1407         /* To prevent confusion, do not again issue a gmx_fatal here since we already
1408          * get the error message from mdrun itself */
1409         sprintf(msg,
1410                 "Cannot run the first benchmark simulation! Please check the error message of\n"
1411                 "mdrun for the source of the problem. Did you provide a command line\n"
1412                 "argument that neither gmx tune_pme nor mdrun understands? If you're\n"
1413                 "sure your command line should work, you can bypass this check with \n"
1414                 "gmx tune_pme -nocheck. The failing command was:\n"
1415                 "\n%s\n\n",
1416                 command);
1417
1418         fprintf(stderr, "%s", msg);
1419         sep_line(fp);
1420         fprintf(fp, "%s", msg);
1421
1422         exit(ret);
1423     }
1424     fprintf(stdout, "Benchmarks can be executed!\n");
1425
1426     /* Clean up the benchmark output files we just created */
1427     fprintf(stdout, "Cleaning up ...\n");
1428     remove_if_exists(opt2fn("-bc", nfile, fnm));
1429     remove_if_exists(opt2fn("-be", nfile, fnm));
1430     remove_if_exists(opt2fn("-bcpo", nfile, fnm));
1431     remove_if_exists(opt2fn("-bg", nfile, fnm));
1432     remove_if_exists(opt2fn("-bo", nfile, fnm));
1433     remove_if_exists(opt2fn("-bx", nfile, fnm));
1434
1435     sfree(command);
1436     sfree(msg);
1437 }
1438
1439 static void do_the_tests(FILE*  fp,                    /* General tune_pme output file           */
1440                          char** tpr_names,             /* Filenames of the input files to test   */
1441                          int    maxPMEnodes,           /* Max fraction of nodes to use for PME   */
1442                          int    minPMEnodes,           /* Min fraction of nodes to use for PME   */
1443                          int    npme_fixed,            /* If >= -1, test fixed number of PME
1444                                                         * nodes only                             */
1445                          const char* npmevalues_opt,   /* Which -npme values should be tested    */
1446                          t_perf**    perfdata,         /* Here the performace data is stored     */
1447                          int*        pmeentries,       /* Entries in the nPMEnodes list          */
1448                          int         repeats,          /* Repeat each test this often            */
1449                          int         nnodes,           /* Total number of nodes = nPP + nPME     */
1450                          int         nr_tprs,          /* Total number of tpr files to test      */
1451                          gmx_bool    bThreads,         /* Threads or MPI?                        */
1452                          char*       cmd_mpirun,       /* mpirun command string                  */
1453                          char*       cmd_np,           /* "-np", "-n", whatever mpirun needs     */
1454                          char*       cmd_mdrun,        /* mdrun command string                   */
1455                          char*       cmd_args_bench,   /* arguments for mdrun in a string        */
1456                          const t_filenm* fnm,          /* List of filenames from command line    */
1457                          int             nfile,        /* Number of files specified on the cmdl. */
1458                          int             presteps,     /* DLB equilibration steps, is checked    */
1459                          int64_t         cpt_steps,    /* Time step counter in the checkpoint    */
1460                          gmx_bool        bCheck,       /* Check whether benchmark mdrun works    */
1461                          const char*     eligible_gpu_ids) /* GPU IDs for
1462                                                         * constructing mdrun command lines */
1463 {
1464     int      i, nr, k, ret, count = 0, totaltests;
1465     int*     nPMEnodes = nullptr;
1466     t_perf*  pd        = nullptr;
1467     int      cmdline_length;
1468     char *   command, *cmd_stub;
1469     char     buf[STRLEN];
1470     gmx_bool bResetProblem = FALSE;
1471     gmx_bool bFirst        = TRUE;
1472
1473     /* This string array corresponds to the eParselog enum type at the start
1474      * of this file */
1475     const char* ParseLog[] = { "OK.",
1476                                "Logfile not found!",
1477                                "No timings, logfile truncated?",
1478                                "Run was terminated.",
1479                                "Counters were not reset properly.",
1480                                "No DD grid found for these settings.",
1481                                "TPX version conflict!",
1482                                "mdrun was not started in parallel!",
1483                                "Number of PP ranks has a prime factor that is too large.",
1484                                "The number of PP ranks did not suit the number of GPUs.",
1485                                "Some GPUs were not detected or are incompatible.",
1486                                "An error occurred." };
1487     char        str_PME_f_load[13];
1488
1489
1490     /* Allocate space for the mdrun command line. 100 extra characters should
1491        be more than enough for the -npme etcetera arguments */
1492     cmdline_length = std::strlen(cmd_mpirun) + std::strlen(cmd_np) + std::strlen(cmd_mdrun)
1493                      + std::strlen(cmd_args_bench) + std::strlen(tpr_names[0]) + 100;
1494     snew(command, cmdline_length);
1495     snew(cmd_stub, cmdline_length);
1496
1497     /* Construct the part of the command line that stays the same for all tests: */
1498     if (bThreads)
1499     {
1500         sprintf(cmd_stub, "%s%s", cmd_mdrun, cmd_np);
1501     }
1502     else
1503     {
1504         sprintf(cmd_stub, "%s%s%s ", cmd_mpirun, cmd_np, cmd_mdrun);
1505     }
1506
1507     /* Create a list of numbers of PME nodes to test */
1508     if (npme_fixed < -1)
1509     {
1510         make_npme_list(npmevalues_opt, pmeentries, &nPMEnodes, nnodes, minPMEnodes, maxPMEnodes);
1511     }
1512     else
1513     {
1514         *pmeentries = 1;
1515         snew(nPMEnodes, 1);
1516         nPMEnodes[0] = npme_fixed;
1517         fprintf(stderr, "Will use a fixed number of %d PME-only ranks.\n", nPMEnodes[0]);
1518     }
1519
1520     if (0 == repeats)
1521     {
1522         fprintf(fp, "\nNo benchmarks done since number of repeats (-r) is 0.\n");
1523         gmx_ffclose(fp);
1524         finalize(opt2fn("-p", nfile, fnm));
1525         exit(0);
1526     }
1527
1528     /* Allocate one dataset for each tpr input file: */
1529     init_perfdata(perfdata, nr_tprs, *pmeentries, repeats);
1530
1531     /*****************************************/
1532     /* Main loop over all tpr files to test: */
1533     /*****************************************/
1534     totaltests = nr_tprs * (*pmeentries) * repeats;
1535     for (k = 0; k < nr_tprs; k++)
1536     {
1537         fprintf(fp, "\nIndividual timings for input file %d (%s):\n", k, tpr_names[k]);
1538         fprintf(fp, "PME ranks      Gcycles       ns/day        PME/f    Remark\n");
1539         /* Loop over various numbers of PME nodes: */
1540         for (i = 0; i < *pmeentries; i++)
1541         {
1542             pd = &perfdata[k][i];
1543
1544             auto cmd_gpu_ids = make_gpu_id_command_line(eligible_gpu_ids);
1545
1546             /* Loop over the repeats for each scenario: */
1547             for (nr = 0; nr < repeats; nr++)
1548             {
1549                 pd->nPMEnodes = nPMEnodes[i];
1550
1551                 /* Add -npme and -s to the command line and save it. Note that
1552                  * the -passall (if set) options requires cmd_args_bench to be
1553                  * at the end of the command line string */
1554                 snew(pd->mdrun_cmd_line, cmdline_length);
1555                 sprintf(pd->mdrun_cmd_line,
1556                         "%s-npme %d -s %s %s %s",
1557                         cmd_stub,
1558                         pd->nPMEnodes,
1559                         tpr_names[k],
1560                         cmd_args_bench,
1561                         cmd_gpu_ids.c_str());
1562
1563                 /* To prevent that all benchmarks fail due to a show-stopper argument
1564                  * on the mdrun command line, we make a quick check first.
1565                  * This check can be turned off in cases where the automatically chosen
1566                  * number of PME-only ranks leads to a number of PP ranks for which no
1567                  * decomposition can be found (e.g. for large prime numbers) */
1568                 if (bFirst && bCheck)
1569                 {
1570                     /* TODO This check is really for a functional
1571                      * .tpr, and if we need it, it should take place
1572                      * for every .tpr, and the logic for it should be
1573                      * immediately inside the loop over k, not in
1574                      * this inner loop. */
1575                     char* temporary_cmd_line;
1576
1577                     snew(temporary_cmd_line, cmdline_length);
1578                     /* TODO -npme 0 is more likely to succeed at low
1579                        parallelism than the default of -npme -1, but
1580                        is more likely to fail at the scaling limit
1581                        when the PP domains may be too small. "mpirun
1582                        -np 1 mdrun" is probably a reasonable thing to
1583                        do for this check, but it'll be easier to
1584                        implement that after some refactoring of how
1585                        the number of MPI ranks is managed. */
1586                     sprintf(temporary_cmd_line, "%s-npme 0 -nb cpu -s %s %s", cmd_stub, tpr_names[k], cmd_args_bench);
1587                     make_sure_it_runs(temporary_cmd_line, cmdline_length, fp, fnm, nfile);
1588                 }
1589                 bFirst = FALSE;
1590
1591                 /* Do a benchmark simulation: */
1592                 if (repeats > 1)
1593                 {
1594                     sprintf(buf, ", pass %d/%d", nr + 1, repeats);
1595                 }
1596                 else
1597                 {
1598                     buf[0] = '\0';
1599                 }
1600                 fprintf(stdout,
1601                         "\n=== Progress %2.0f%%, tpr %d/%d, run %d/%d%s:\n",
1602                         (100.0 * count) / totaltests,
1603                         k + 1,
1604                         nr_tprs,
1605                         i + 1,
1606                         *pmeentries,
1607                         buf);
1608                 make_backup(opt2fn("-err", nfile, fnm));
1609                 sprintf(command, "%s 1> /dev/null 2>%s", pd->mdrun_cmd_line, opt2fn("-err", nfile, fnm));
1610                 fprintf(stdout, "%s\n", pd->mdrun_cmd_line);
1611                 gmx_system_call(command);
1612
1613                 /* Collect the performance data from the log file; also check stderr
1614                  * for fatal errors */
1615                 ret = parse_logfile(
1616                         opt2fn("-bg", nfile, fnm), opt2fn("-err", nfile, fnm), pd, nr, presteps, cpt_steps, nnodes);
1617                 if ((presteps > 0) && (ret == eParselogResetProblem))
1618                 {
1619                     bResetProblem = TRUE;
1620                 }
1621
1622                 if (-1 == pd->nPMEnodes)
1623                 {
1624                     sprintf(buf, "(%3d)", pd->guessPME);
1625                 }
1626                 else
1627                 {
1628                     sprintf(buf, "     ");
1629                 }
1630
1631                 /* Nicer output */
1632                 if (pd->PME_f_load[nr] > 0.0)
1633                 {
1634                     sprintf(str_PME_f_load, "%12.3f", pd->PME_f_load[nr]);
1635                 }
1636                 else
1637                 {
1638                     sprintf(str_PME_f_load, "%s", "         -  ");
1639                 }
1640
1641                 /* Write the data we got to disk */
1642                 fprintf(fp,
1643                         "%4d%s %12.3f %12.3f %s    %s",
1644                         pd->nPMEnodes,
1645                         buf,
1646                         pd->Gcycles[nr],
1647                         pd->ns_per_day[nr],
1648                         str_PME_f_load,
1649                         ParseLog[ret]);
1650                 if (!(ret == eParselogOK || ret == eParselogNoDDGrid || ret == eParselogNotFound))
1651                 {
1652                     fprintf(fp, " Check %s file for problems.", ret == eParselogFatal ? "err" : "log");
1653                 }
1654                 fprintf(fp, "\n");
1655                 fflush(fp);
1656                 count++;
1657
1658                 /* Do some cleaning up and delete the files we do not need any more */
1659                 cleanup(fnm, nfile, k, nnodes, pd->nPMEnodes, nr, ret == eParselogFatal);
1660
1661                 /* If the first run with this number of processors already failed, do not try again: */
1662                 if (pd->Gcycles[0] <= 0.0 && repeats > 1)
1663                 {
1664                     fprintf(stdout,
1665                             "Skipping remaining passes of unsuccessful setting, see log file for "
1666                             "details.\n");
1667                     count += repeats - (nr + 1);
1668                     break;
1669                 }
1670             } /* end of repeats loop */
1671         }     /* end of -npme loop */
1672     }         /* end of tpr file loop */
1673
1674     if (bResetProblem)
1675     {
1676         sep_line(fp);
1677         fprintf(fp, "WARNING: The cycle and time step counters could not be reset properly. ");
1678         sep_line(fp);
1679     }
1680     sfree(command);
1681     sfree(cmd_stub);
1682 }
1683
1684
1685 static void check_input(int             nnodes,
1686                         int             repeats,
1687                         int*            ntprs,
1688                         real*           rmin,
1689                         real            rcoulomb,
1690                         real*           rmax,
1691                         real            maxPMEfraction,
1692                         real            minPMEfraction,
1693                         int             npme_fixed,
1694                         int64_t         bench_nsteps,
1695                         const t_filenm* fnm,
1696                         int             nfile,
1697                         int             sim_part,
1698                         int             presteps,
1699                         int             npargs,
1700                         t_pargs*        pa)
1701 {
1702     int old;
1703
1704
1705     /* Make sure the input file exists */
1706     if (!gmx_fexist(opt2fn("-s", nfile, fnm)))
1707     {
1708         gmx_fatal(FARGS, "File %s not found.", opt2fn("-s", nfile, fnm));
1709     }
1710
1711     /* Make sure that the checkpoint file is not overwritten during benchmarking */
1712     if ((0 == std::strcmp(opt2fn("-cpi", nfile, fnm), opt2fn("-bcpo", nfile, fnm))) && (sim_part > 1))
1713     {
1714         gmx_fatal(FARGS,
1715                   "Checkpoint input (-cpi) and benchmark checkpoint output (-bcpo) files must not "
1716                   "be identical.\n"
1717                   "The checkpoint input file must not be overwritten during the benchmarks.\n");
1718     }
1719
1720     /* Make sure that repeats is >= 0 (if == 0, only write tpr files) */
1721     if (repeats < 0)
1722     {
1723         gmx_fatal(FARGS, "Number of repeats < 0!");
1724     }
1725
1726     /* Check number of nodes */
1727     if (nnodes < 1)
1728     {
1729         gmx_fatal(FARGS, "Number of ranks/threads must be a positive integer.");
1730     }
1731
1732     /* Automatically choose -ntpr if not set */
1733     if (*ntprs < 1)
1734     {
1735         if (nnodes < 16)
1736         {
1737             *ntprs = 1;
1738         }
1739         else
1740         {
1741             *ntprs = 3;
1742             /* Set a reasonable scaling factor for rcoulomb */
1743             if (*rmax <= 0)
1744             {
1745                 *rmax = rcoulomb * 1.2;
1746             }
1747         }
1748         fprintf(stderr, "Will test %d tpr file%s.\n", *ntprs, *ntprs == 1 ? "" : "s");
1749     }
1750     else
1751     {
1752         if (1 == *ntprs)
1753         {
1754             fprintf(stderr,
1755                     "Note: Choose ntpr>1 to shift PME load between real and reciprocal space.\n");
1756         }
1757     }
1758
1759     /* Make shure that rmin <= rcoulomb <= rmax */
1760     if (*rmin <= 0)
1761     {
1762         *rmin = rcoulomb;
1763     }
1764     if (*rmax <= 0)
1765     {
1766         *rmax = rcoulomb;
1767     }
1768     if (!(*rmin <= *rmax))
1769     {
1770         gmx_fatal(FARGS,
1771                   "Please choose the Coulomb radii such that rmin <= rmax.\n"
1772                   "rmin = %g, rmax = %g, actual rcoul from .tpr file = %g\n",
1773                   *rmin,
1774                   *rmax,
1775                   rcoulomb);
1776     }
1777     /* Add test scenarios if rmin or rmax were set */
1778     if (*ntprs <= 2)
1779     {
1780         if (!gmx_within_tol(*rmin, rcoulomb, GMX_REAL_EPS) && (*ntprs == 1))
1781         {
1782             (*ntprs)++;
1783             fprintf(stderr, "NOTE: Setting -rmin to %g changed -ntpr to %d\n", *rmin, *ntprs);
1784         }
1785         if (!gmx_within_tol(*rmax, rcoulomb, GMX_REAL_EPS) && (*ntprs == 1))
1786         {
1787             (*ntprs)++;
1788             fprintf(stderr, "NOTE: Setting -rmax to %g changed -ntpr to %d\n", *rmax, *ntprs);
1789         }
1790     }
1791     old = *ntprs;
1792     /* If one of rmin, rmax is set, we need 2 tpr files at minimum */
1793     if (!gmx_within_tol(*rmax, rcoulomb, GMX_REAL_EPS) || !gmx_within_tol(*rmin, rcoulomb, GMX_REAL_EPS))
1794     {
1795         *ntprs = std::max(*ntprs, 2);
1796     }
1797
1798     /* If both rmin, rmax are set, we need 3 tpr files at minimum */
1799     if (!gmx_within_tol(*rmax, rcoulomb, GMX_REAL_EPS) && !gmx_within_tol(*rmin, rcoulomb, GMX_REAL_EPS))
1800     {
1801         *ntprs = std::max(*ntprs, 3);
1802     }
1803
1804     if (old != *ntprs)
1805     {
1806         fprintf(stderr, "NOTE: Your rmin, rmax setting changed -ntpr to %d\n", *ntprs);
1807     }
1808
1809     if (*ntprs > 1)
1810     {
1811         if (gmx_within_tol(*rmin, rcoulomb, GMX_REAL_EPS)
1812             && gmx_within_tol(rcoulomb, *rmax, GMX_REAL_EPS)) /* We have just a single rc */
1813         {
1814             fprintf(stderr,
1815                     "WARNING: Resetting -ntpr to 1 since no Coulomb radius scaling is requested.\n"
1816                     "Please set rmin < rmax to test Coulomb radii in the [rmin, rmax] interval\n"
1817                     "with correspondingly adjusted PME grid settings\n");
1818             *ntprs = 1;
1819         }
1820     }
1821
1822     /* Check whether max and min fraction are within required values */
1823     if (maxPMEfraction > 0.5 || maxPMEfraction < 0)
1824     {
1825         gmx_fatal(FARGS, "-max must be between 0 and 0.5");
1826     }
1827     if (minPMEfraction > 0.5 || minPMEfraction < 0)
1828     {
1829         gmx_fatal(FARGS, "-min must be between 0 and 0.5");
1830     }
1831     if (maxPMEfraction < minPMEfraction)
1832     {
1833         gmx_fatal(FARGS, "-max must be larger or equal to -min");
1834     }
1835
1836     /* Check whether the number of steps - if it was set - has a reasonable value */
1837     if (bench_nsteps < 0)
1838     {
1839         gmx_fatal(FARGS, "Number of steps must be positive.");
1840     }
1841
1842     if (bench_nsteps > 10000 || bench_nsteps < 100)
1843     {
1844         fprintf(stderr, "WARNING: steps=");
1845         fprintf(stderr, "%" PRId64, bench_nsteps);
1846         fprintf(stderr,
1847                 ". Are you sure you want to perform so %s steps for each benchmark?\n",
1848                 (bench_nsteps < 100) ? "few" : "many");
1849     }
1850
1851     if (presteps < 0)
1852     {
1853         gmx_fatal(FARGS, "Cannot have a negative number of presteps.\n");
1854     }
1855
1856     /* Check for rcoulomb scaling if more than one .tpr file is tested */
1857     if (*ntprs > 1)
1858     {
1859         if (*rmin / rcoulomb < 0.75 || *rmax / rcoulomb > 1.25)
1860         {
1861             fprintf(stderr,
1862                     "WARNING: Applying extreme scaling factor. I hope you know what you are "
1863                     "doing.\n");
1864         }
1865     }
1866
1867     /* If a fixed number of PME nodes is set we do rcoulomb and PME gird tuning
1868      * only. We need to check whether the requested number of PME-only nodes
1869      * makes sense. */
1870     if (npme_fixed > -1)
1871     {
1872         /* No more than 50% of all nodes can be assigned as PME-only nodes. */
1873         if (2 * npme_fixed > nnodes)
1874         {
1875             gmx_fatal(FARGS,
1876                       "Cannot have more than %d PME-only ranks for a total of %d ranks (you chose "
1877                       "%d).\n",
1878                       nnodes / 2,
1879                       nnodes,
1880                       npme_fixed);
1881         }
1882         if ((npme_fixed > 0) && (5 * npme_fixed < nnodes))
1883         {
1884             fprintf(stderr,
1885                     "WARNING: Only %g percent of the ranks are assigned as PME-only ranks.\n",
1886                     (100.0 * npme_fixed) / nnodes);
1887         }
1888         if (opt2parg_bSet("-min", npargs, pa) || opt2parg_bSet("-max", npargs, pa))
1889         {
1890             fprintf(stderr,
1891                     "NOTE: The -min, -max, and -npme options have no effect when a\n"
1892                     "      fixed number of PME-only ranks is requested with -fix.\n");
1893         }
1894     }
1895 }
1896
1897
1898 /* Returns TRUE when "opt" is needed at launch time */
1899 static gmx_bool is_launch_file(char* opt, gmx_bool bSet)
1900 {
1901     if (0 == std::strncmp(opt, "-swap", 5))
1902     {
1903         return bSet;
1904     }
1905
1906     /* Apart from the input .tpr and the output log files we need all options that
1907      * were set on the command line and that do not start with -b */
1908     if (0 == std::strncmp(opt, "-b", 2) || 0 == std::strncmp(opt, "-s", 2)
1909         || 0 == std::strncmp(opt, "-err", 4) || 0 == std::strncmp(opt, "-p", 2))
1910     {
1911         return FALSE;
1912     }
1913
1914     return bSet;
1915 }
1916
1917
1918 /* Returns TRUE when "opt" defines a file which is needed for the benchmarks runs */
1919 static gmx_bool is_bench_file(char* opt, gmx_bool bSet, gmx_bool bOptional, gmx_bool bIsOutput)
1920 {
1921     /* Apart from the input .tpr, all files starting with "-b" are for
1922      * _b_enchmark files exclusively */
1923     if (0 == std::strncmp(opt, "-s", 2))
1924     {
1925         return FALSE;
1926     }
1927
1928     if (0 == std::strncmp(opt, "-b", 2) || 0 == std::strncmp(opt, "-s", 2))
1929     {
1930         return !bOptional || bSet;
1931     }
1932     else
1933     {
1934         if (bIsOutput)
1935         {
1936             return FALSE;
1937         }
1938         else
1939         {
1940             return bSet; /* These are additional input files like -cpi -ei */
1941         }
1942     }
1943 }
1944
1945
1946 /* Adds 'buf' to 'str' */
1947 static void add_to_string(char** str, const char* buf)
1948 {
1949     int len;
1950
1951
1952     len = std::strlen(*str) + std::strlen(buf) + 1;
1953     srenew(*str, len);
1954     std::strcat(*str, buf);
1955 }
1956
1957
1958 /* Create the command line for the benchmark as well as for the real run */
1959 static void
1960 create_command_line_snippets(gmx_bool bAppendFiles,
1961                              gmx_bool bKeepAndNumCPT,
1962                              gmx_bool bResetHWay,
1963                              int      presteps,
1964                              int      nfile,
1965                              t_filenm fnm[],
1966                              char* cmd_args_bench[], /* command line arguments for benchmark runs */
1967                              char* cmd_args_launch[], /* command line arguments for simulation run */
1968                              char  extra_args[],      /* Add this to the end of the command line */
1969                              char* deffnm)            /* Default file names, or NULL if not set */
1970 {
1971     int         i;
1972     char*       opt;
1973     const char* name;
1974     char        strbuf[STRLEN];
1975
1976
1977     /* strlen needs at least '\0' as a string: */
1978     snew(*cmd_args_bench, 1);
1979     snew(*cmd_args_launch, 1);
1980     *cmd_args_launch[0] = '\0';
1981     *cmd_args_bench[0]  = '\0';
1982
1983
1984     /*******************************************/
1985     /* 1. Process other command line arguments */
1986     /*******************************************/
1987     if (presteps > 0)
1988     {
1989         /* Add equilibration steps to benchmark options */
1990         sprintf(strbuf, "-resetstep %d ", presteps);
1991         add_to_string(cmd_args_bench, strbuf);
1992     }
1993     /* These switches take effect only at launch time */
1994     if (deffnm)
1995     {
1996         sprintf(strbuf, "-deffnm %s ", deffnm);
1997         add_to_string(cmd_args_launch, strbuf);
1998     }
1999     if (!bAppendFiles)
2000     {
2001         add_to_string(cmd_args_launch, "-noappend ");
2002     }
2003     if (bKeepAndNumCPT)
2004     {
2005         add_to_string(cmd_args_launch, "-cpnum ");
2006     }
2007     if (bResetHWay)
2008     {
2009         add_to_string(cmd_args_launch, "-resethway ");
2010     }
2011
2012     /********************/
2013     /* 2. Process files */
2014     /********************/
2015     for (i = 0; i < nfile; i++)
2016     {
2017         opt  = const_cast<char*>(fnm[i].opt);
2018         name = opt2fn(opt, nfile, fnm);
2019
2020         /* Strbuf contains the options, now let's sort out where we need that */
2021         sprintf(strbuf, "%s %s ", opt, name);
2022
2023         if (is_bench_file(opt, opt2bSet(opt, nfile, fnm), is_optional(&fnm[i]), is_output(&fnm[i])))
2024         {
2025             /* All options starting with -b* need the 'b' removed,
2026              * therefore overwrite strbuf */
2027             if (0 == std::strncmp(opt, "-b", 2))
2028             {
2029                 sprintf(strbuf, "-%s %s ", &opt[2], name);
2030             }
2031
2032             add_to_string(cmd_args_bench, strbuf);
2033         }
2034
2035         if (is_launch_file(opt, opt2bSet(opt, nfile, fnm)))
2036         {
2037             add_to_string(cmd_args_launch, strbuf);
2038         }
2039     }
2040
2041     add_to_string(cmd_args_bench, extra_args);
2042     add_to_string(cmd_args_launch, extra_args);
2043 }
2044
2045
2046 /* Set option opt */
2047 static void setopt(const char* opt, int nfile, t_filenm fnm[])
2048 {
2049     int i;
2050
2051     for (i = 0; (i < nfile); i++)
2052     {
2053         if (std::strcmp(opt, fnm[i].opt) == 0)
2054         {
2055             fnm[i].flag |= ffSET;
2056         }
2057     }
2058 }
2059
2060
2061 /* This routine inspects the tpr file and ...
2062  * 1. checks for output files that get triggered by a tpr option. These output
2063  *    files are marked as 'set' to allow for a proper cleanup after each
2064  *    tuning run.
2065  * 2. returns the PME:PP load ratio
2066  * 3. returns rcoulomb from the tpr */
2067 static float inspect_tpr(int nfile, t_filenm fnm[], real* rcoulomb)
2068 {
2069     gmx_bool   bTpi;  /* Is test particle insertion requested?          */
2070     gmx_bool   bFree; /* Is a free energy simulation requested?         */
2071     gmx_bool   bNM;   /* Is a normal mode analysis requested?           */
2072     gmx_bool   bSwap; /* Is water/ion position swapping requested?      */
2073     t_state    state;
2074     gmx_mtop_t mtop;
2075
2076
2077     /* Check tpr file for options that trigger extra output files */
2078     t_inputrec  irInstance;
2079     t_inputrec* ir = &irInstance;
2080     read_tpx_state(opt2fn("-s", nfile, fnm), ir, &state, &mtop);
2081     bFree = (FreeEnergyPerturbationType::No != ir->efep);
2082     bNM   = (IntegrationAlgorithm::NM == ir->eI);
2083     bSwap = (SwapType::No != ir->eSwapCoords);
2084     bTpi  = EI_TPI(ir->eI);
2085
2086     /* Set these output files on the tuning command-line */
2087     if (ir->bPull)
2088     {
2089         setopt("-pf", nfile, fnm);
2090         setopt("-px", nfile, fnm);
2091     }
2092     if (bFree)
2093     {
2094         setopt("-dhdl", nfile, fnm);
2095     }
2096     if (bTpi)
2097     {
2098         setopt("-tpi", nfile, fnm);
2099         setopt("-tpid", nfile, fnm);
2100     }
2101     if (bNM)
2102     {
2103         setopt("-mtx", nfile, fnm);
2104     }
2105     if (bSwap)
2106     {
2107         setopt("-swap", nfile, fnm);
2108     }
2109
2110     *rcoulomb = ir->rcoulomb;
2111
2112     /* Return the estimate for the number of PME nodes */
2113     float npme = pme_load_estimate(mtop, *ir, state.box);
2114     return npme;
2115 }
2116
2117
2118 static void couple_files_options(int nfile, t_filenm fnm[])
2119 {
2120     int      i;
2121     gmx_bool bSet, bBench;
2122     char*    opt;
2123     char     buf[20];
2124
2125
2126     for (i = 0; i < nfile; i++)
2127     {
2128         opt    = const_cast<char*>(fnm[i].opt);
2129         bSet   = ((fnm[i].flag & ffSET) != 0);
2130         bBench = (0 == std::strncmp(opt, "-b", 2));
2131
2132         /* Check optional files */
2133         /* If e.g. -eo is set, then -beo also needs to be set */
2134         if (is_optional(&fnm[i]) && bSet && !bBench)
2135         {
2136             sprintf(buf, "-b%s", &opt[1]);
2137             setopt(buf, nfile, fnm);
2138         }
2139         /* If -beo is set, then -eo also needs to be! */
2140         if (is_optional(&fnm[i]) && bSet && bBench)
2141         {
2142             sprintf(buf, "-%s", &opt[2]);
2143             setopt(buf, nfile, fnm);
2144         }
2145     }
2146 }
2147
2148
2149 #define BENCHSTEPS (1000)
2150
2151 int gmx_tune_pme(int argc, char* argv[])
2152 {
2153     const char* desc[] = {
2154         "For a given number [TT]-np[tt] or [TT]-ntmpi[tt] of ranks, [THISMODULE] systematically",
2155         "times [gmx-mdrun] with various numbers of PME-only ranks and determines",
2156         "which setting is fastest. It will also test whether performance can",
2157         "be enhanced by shifting load from the reciprocal to the real space",
2158         "part of the Ewald sum. ",
2159         "Simply pass your [REF].tpr[ref] file to [THISMODULE] together with other options",
2160         "for [gmx-mdrun] as needed.[PAR]",
2161         "[THISMODULE] needs to call [gmx-mdrun] and so requires that you",
2162         "specify how to call mdrun with the argument to the [TT]-mdrun[tt]",
2163         "parameter. Depending how you have built GROMACS, values such as",
2164         "'gmx mdrun', 'gmx_d mdrun', or 'gmx_mpi mdrun' might be needed.[PAR]",
2165         "The program that runs MPI programs can be set in the environment variable",
2166         "MPIRUN (defaults to 'mpirun'). Note that for certain MPI frameworks,",
2167         "you need to provide a machine- or hostfile. This can also be passed",
2168         "via the MPIRUN variable, e.g.[PAR]",
2169         "[TT]export MPIRUN=\"/usr/local/mpirun -machinefile hosts\"[tt]",
2170         "Note that in such cases it is normally necessary to compile",
2171         "and/or run [THISMODULE] without MPI support, so that it can call",
2172         "the MPIRUN program.[PAR]",
2173         "Before doing the actual benchmark runs, [THISMODULE] will do a quick",
2174         "check whether [gmx-mdrun] works as expected with the provided parallel settings",
2175         "if the [TT]-check[tt] option is activated (the default).",
2176         "Please call [THISMODULE] with the normal options you would pass to",
2177         "[gmx-mdrun] and add [TT]-np[tt] for the number of ranks to perform the",
2178         "tests on, or [TT]-ntmpi[tt] for the number of threads. You can also add [TT]-r[tt]",
2179         "to repeat each test several times to get better statistics. [PAR]",
2180         "[THISMODULE] can test various real space / reciprocal space workloads",
2181         "for you. With [TT]-ntpr[tt] you control how many extra [REF].tpr[ref] files will be",
2182         "written with enlarged cutoffs and smaller Fourier grids respectively.",
2183         "Typically, the first test (number 0) will be with the settings from the input",
2184         "[REF].tpr[ref] file; the last test (number [TT]ntpr[tt]) will have the Coulomb cutoff",
2185         "specified by [TT]-rmax[tt] with a somewhat smaller PME grid at the same time. ",
2186         "In this last test, the Fourier spacing is multiplied with [TT]rmax[tt]/rcoulomb. ",
2187         "The remaining [REF].tpr[ref] files will have equally-spaced Coulomb radii (and Fourier ",
2188         "spacings) between these extremes. [BB]Note[bb] that you can set [TT]-ntpr[tt] to 1",
2189         "if you just seek the optimal number of PME-only ranks; in that case",
2190         "your input [REF].tpr[ref] file will remain unchanged.[PAR]",
2191         "For the benchmark runs, the default of 1000 time steps should suffice for most",
2192         "MD systems. The dynamic load balancing needs about 100 time steps",
2193         "to adapt to local load imbalances, therefore the time step counters",
2194         "are by default reset after 100 steps. For large systems (>1M atoms), as well as ",
2195         "for a higher accuracy of the measurements, you should set [TT]-resetstep[tt] to a higher ",
2196         "value.",
2197         "From the 'DD' load imbalance entries in the md.log output file you",
2198         "can tell after how many steps the load is sufficiently balanced. Example call:[PAR]",
2199         "[TT]gmx tune_pme -np 64 -s protein.tpr -launch[tt][PAR]",
2200         "After calling [gmx-mdrun] several times, detailed performance information",
2201         "is available in the output file [TT]perf.out[tt].",
2202         "[BB]Note[bb] that during the benchmarks, a couple of temporary files are written",
2203         "(options [TT]-b*[tt]), these will be automatically deleted after each test.[PAR]",
2204         "If you want the simulation to be started automatically with the",
2205         "optimized parameters, use the command line option [TT]-launch[tt].[PAR]",
2206         "Basic support for GPU-enabled [TT]mdrun[tt] exists. Give a string containing the IDs",
2207         "of the GPUs that you wish to use in the optimization in the [TT]-gpu_id[tt]",
2208         "command-line argument. This works exactly like [TT]mdrun -gpu_id[tt], does not imply a ",
2209         "mapping,",
2210         "and merely declares the eligible set of GPU devices. [TT]gmx-tune_pme[tt] will construct ",
2211         "calls to",
2212         "mdrun that use this set appropriately. [TT]gmx-tune_pme[tt] does not support",
2213         "[TT]-gputasks[tt].[PAR]",
2214     };
2215
2216     int   nnodes         = 1;
2217     int   repeats        = 2;
2218     int   pmeentries     = 0; /* How many values for -npme do we actually test for each tpr file */
2219     real  maxPMEfraction = 0.50;
2220     real  minPMEfraction = 0.25;
2221     int   maxPMEnodes, minPMEnodes;
2222     float guessPMEratio; /* guessed PME:PP ratio based on the tpr file */
2223     float guessPMEnodes;
2224     int   npme_fixed = -2; /* If >= -1, use only this number
2225                             * of PME-only nodes                */
2226     int         ntprs = 0;
2227     real        rmin = 0.0, rmax = 0.0; /* min and max value for rcoulomb if scaling is requested */
2228     real        rcoulomb       = -1.0;  /* Coulomb radius as set in .tpr file */
2229     gmx_bool    bScaleRvdw     = TRUE;
2230     int64_t     bench_nsteps   = BENCHSTEPS;
2231     int64_t     new_sim_nsteps = -1;   /* -1 indicates: not set by the user */
2232     int64_t     cpt_steps      = 0;    /* Step counter in .cpt input file   */
2233     int         presteps       = 1500; /* Do a full cycle reset after presteps steps */
2234     gmx_bool    bOverwrite     = FALSE, bKeepTPR;
2235     gmx_bool    bLaunch        = FALSE;
2236     char*       ExtraArgs      = nullptr;
2237     char**      tpr_names      = nullptr;
2238     const char* simulation_tpr = nullptr;
2239     char*       deffnm         = nullptr;
2240     int         best_npme, best_tpr;
2241     int         sim_part = 1; /* For benchmarks with checkpoint files */
2242     char        bbuf[STRLEN];
2243
2244
2245     /* Default program names if nothing else is found */
2246     char *cmd_mpirun = nullptr, *cmd_mdrun = nullptr;
2247     char *cmd_args_bench, *cmd_args_launch;
2248     char* cmd_np = nullptr;
2249
2250     /* IDs of GPUs that are eligible for computation */
2251     char* eligible_gpu_ids = nullptr;
2252
2253     t_perf**     perfdata = nullptr;
2254     t_inputinfo* info;
2255     int          i;
2256     FILE*        fp;
2257
2258     /* Print out how long the tuning took */
2259     double seconds;
2260
2261     static t_filenm fnm[] = { /* tune_pme */
2262                               { efOUT, "-p", "perf", ffWRITE },
2263                               { efLOG, "-err", "bencherr", ffWRITE },
2264                               { efTPR, "-so", "tuned", ffWRITE },
2265                               /* mdrun: */
2266                               { efTPR, "-s", nullptr, ffREAD },
2267                               { efTRN, "-o", nullptr, ffWRITE },
2268                               { efCOMPRESSED, "-x", nullptr, ffOPTWR },
2269                               { efCPT, "-cpi", nullptr, ffOPTRD },
2270                               { efCPT, "-cpo", nullptr, ffOPTWR },
2271                               { efSTO, "-c", "confout", ffWRITE },
2272                               { efEDR, "-e", "ener", ffWRITE },
2273                               { efLOG, "-g", "md", ffWRITE },
2274                               { efXVG, "-dhdl", "dhdl", ffOPTWR },
2275                               { efXVG, "-field", "field", ffOPTWR },
2276                               { efXVG, "-table", "table", ffOPTRD },
2277                               { efXVG, "-tablep", "tablep", ffOPTRD },
2278                               { efXVG, "-tableb", "table", ffOPTRD },
2279                               { efTRX, "-rerun", "rerun", ffOPTRD },
2280                               { efXVG, "-tpi", "tpi", ffOPTWR },
2281                               { efXVG, "-tpid", "tpidist", ffOPTWR },
2282                               { efEDI, "-ei", "sam", ffOPTRD },
2283                               { efXVG, "-eo", "edsam", ffOPTWR },
2284                               { efXVG, "-px", "pullx", ffOPTWR },
2285                               { efXVG, "-pf", "pullf", ffOPTWR },
2286                               { efXVG, "-ro", "rotation", ffOPTWR },
2287                               { efLOG, "-ra", "rotangles", ffOPTWR },
2288                               { efLOG, "-rs", "rotslabs", ffOPTWR },
2289                               { efLOG, "-rt", "rottorque", ffOPTWR },
2290                               { efMTX, "-mtx", "nm", ffOPTWR },
2291                               { efXVG, "-swap", "swapions", ffOPTWR },
2292                               /* Output files that are deleted after each benchmark run */
2293                               { efTRN, "-bo", "bench", ffWRITE },
2294                               { efXTC, "-bx", "bench", ffWRITE },
2295                               { efCPT, "-bcpo", "bench", ffWRITE },
2296                               { efSTO, "-bc", "bench", ffWRITE },
2297                               { efEDR, "-be", "bench", ffWRITE },
2298                               { efLOG, "-bg", "bench", ffWRITE },
2299                               { efXVG, "-beo", "benchedo", ffOPTWR },
2300                               { efXVG, "-bdhdl", "benchdhdl", ffOPTWR },
2301                               { efXVG, "-bfield", "benchfld", ffOPTWR },
2302                               { efXVG, "-btpi", "benchtpi", ffOPTWR },
2303                               { efXVG, "-btpid", "benchtpid", ffOPTWR },
2304                               { efXVG, "-bdevout", "benchdev", ffOPTWR },
2305                               { efXVG, "-brunav", "benchrnav", ffOPTWR },
2306                               { efXVG, "-bpx", "benchpx", ffOPTWR },
2307                               { efXVG, "-bpf", "benchpf", ffOPTWR },
2308                               { efXVG, "-bro", "benchrot", ffOPTWR },
2309                               { efLOG, "-bra", "benchrota", ffOPTWR },
2310                               { efLOG, "-brs", "benchrots", ffOPTWR },
2311                               { efLOG, "-brt", "benchrott", ffOPTWR },
2312                               { efMTX, "-bmtx", "benchn", ffOPTWR },
2313                               { efNDX, "-bdn", "bench", ffOPTWR },
2314                               { efXVG, "-bswap", "benchswp", ffOPTWR }
2315     };
2316
2317     gmx_bool bThreads = FALSE;
2318
2319     int nthreads = 1;
2320
2321     const char* procstring[]     = { nullptr, "np", "n", "none", nullptr };
2322     const char* npmevalues_opt[] = { nullptr, "auto", "all", "subset", nullptr };
2323
2324     gmx_bool bAppendFiles          = TRUE;
2325     gmx_bool bKeepAndNumCPT        = FALSE;
2326     gmx_bool bResetCountersHalfWay = FALSE;
2327     gmx_bool bBenchmark            = TRUE;
2328     gmx_bool bCheck                = TRUE;
2329
2330     gmx_output_env_t* oenv = nullptr;
2331
2332     t_pargs pa[] = {
2333         /***********************/
2334         /* tune_pme options: */
2335         /***********************/
2336         { "-mdrun",
2337           FALSE,
2338           etSTR,
2339           { &cmd_mdrun },
2340           "Command line to run a simulation, e.g. 'gmx mdrun' or 'gmx_mpi mdrun'" },
2341         { "-np",
2342           FALSE,
2343           etINT,
2344           { &nnodes },
2345           "Number of ranks to run the tests on (must be > 2 for separate PME ranks)" },
2346         { "-npstring",
2347           FALSE,
2348           etENUM,
2349           { procstring },
2350           "Name of the [TT]$MPIRUN[tt] option that specifies the number of ranks to use ('np', or "
2351           "'n'; use 'none' if there is no such option)" },
2352         { "-ntmpi",
2353           FALSE,
2354           etINT,
2355           { &nthreads },
2356           "Number of MPI-threads to run the tests on (turns MPI & mpirun off)" },
2357         { "-r", FALSE, etINT, { &repeats }, "Repeat each test this often" },
2358         { "-max", FALSE, etREAL, { &maxPMEfraction }, "Max fraction of PME ranks to test with" },
2359         { "-min", FALSE, etREAL, { &minPMEfraction }, "Min fraction of PME ranks to test with" },
2360         { "-npme",
2361           FALSE,
2362           etENUM,
2363           { npmevalues_opt },
2364           "Within -min and -max, benchmark all possible values for [TT]-npme[tt], or just a "
2365           "reasonable subset. "
2366           "Auto neglects -min and -max and chooses reasonable values around a guess for npme "
2367           "derived from the .tpr" },
2368         { "-fix",
2369           FALSE,
2370           etINT,
2371           { &npme_fixed },
2372           "If >= -1, do not vary the number of PME-only ranks, instead use this fixed value and "
2373           "only vary rcoulomb and the PME grid spacing." },
2374         { "-rmax",
2375           FALSE,
2376           etREAL,
2377           { &rmax },
2378           "If >0, maximal rcoulomb for -ntpr>1 (rcoulomb upscaling results in fourier grid "
2379           "downscaling)" },
2380         { "-rmin", FALSE, etREAL, { &rmin }, "If >0, minimal rcoulomb for -ntpr>1" },
2381         { "-scalevdw", FALSE, etBOOL, { &bScaleRvdw }, "Scale rvdw along with rcoulomb" },
2382         { "-ntpr",
2383           FALSE,
2384           etINT,
2385           { &ntprs },
2386           "Number of [REF].tpr[ref] files to benchmark. Create this many files with different "
2387           "rcoulomb scaling factors depending on -rmin and -rmax. "
2388           "If < 1, automatically choose the number of [REF].tpr[ref] files to test" },
2389         { "-steps",
2390           FALSE,
2391           etINT64,
2392           { &bench_nsteps },
2393           "Take timings for this many steps in the benchmark runs" },
2394         { "-resetstep",
2395           FALSE,
2396           etINT,
2397           { &presteps },
2398           "Let dlb equilibrate this many steps before timings are taken (reset cycle counters "
2399           "after this many steps)" },
2400         { "-nsteps",
2401           FALSE,
2402           etINT64,
2403           { &new_sim_nsteps },
2404           "If non-negative, perform this many steps in the real run (overwrites nsteps from "
2405           "[REF].tpr[ref], add [REF].cpt[ref] steps)" },
2406         { "-launch", FALSE, etBOOL, { &bLaunch }, "Launch the real simulation after optimization" },
2407         { "-bench",
2408           FALSE,
2409           etBOOL,
2410           { &bBenchmark },
2411           "Run the benchmarks or just create the input [REF].tpr[ref] files?" },
2412         { "-check",
2413           FALSE,
2414           etBOOL,
2415           { &bCheck },
2416           "Before the benchmark runs, check whether mdrun works in parallel" },
2417         { "-gpu_id",
2418           FALSE,
2419           etSTR,
2420           { &eligible_gpu_ids },
2421           "List of unique GPU device IDs that are eligible for use" },
2422         /******************/
2423         /* mdrun options: */
2424         /******************/
2425         /* We let tune_pme parse and understand these options, because we need to
2426          * prevent that they appear on the mdrun command line for the benchmarks */
2427         { "-append",
2428           FALSE,
2429           etBOOL,
2430           { &bAppendFiles },
2431           "Append to previous output files when continuing from checkpoint instead of adding the "
2432           "simulation part number to all file names (for launch only)" },
2433         { "-cpnum",
2434           FALSE,
2435           etBOOL,
2436           { &bKeepAndNumCPT },
2437           "Keep and number checkpoint files (launch only)" },
2438         { "-deffnm", FALSE, etSTR, { &deffnm }, "Set the default filenames (launch only)" },
2439         { "-resethway",
2440           FALSE,
2441           etBOOL,
2442           { &bResetCountersHalfWay },
2443           "HIDDENReset the cycle counters after half the number of steps or halfway [TT]-maxh[tt] "
2444           "(launch only)" }
2445     };
2446
2447 #define NFILE asize(fnm)
2448
2449     seconds = gmx_gettime();
2450
2451     if (!parse_common_args(
2452                 &argc, argv, PCA_NOEXIT_ON_ARGS, NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
2453     {
2454         return 0;
2455     }
2456
2457     // procstring[0] is used inside two different conditionals further down
2458     GMX_RELEASE_ASSERT(procstring[0] != nullptr, "Options inconsistency; procstring[0] is NULL");
2459
2460     /* Store the remaining unparsed command line entries in a string which
2461      * is then attached to the mdrun command line */
2462     snew(ExtraArgs, 1);
2463     ExtraArgs[0] = '\0';
2464     for (i = 1; i < argc; i++) /* argc will now be 1 if everything was understood */
2465     {
2466         add_to_string(&ExtraArgs, argv[i]);
2467         add_to_string(&ExtraArgs, " ");
2468     }
2469
2470     if (opt2parg_bSet("-ntmpi", asize(pa), pa))
2471     {
2472         bThreads = TRUE;
2473         if (opt2parg_bSet("-npstring", asize(pa), pa))
2474         {
2475             fprintf(stderr, "WARNING: -npstring has no effect when using threads.\n");
2476         }
2477
2478         if (nnodes > 1)
2479         {
2480             gmx_fatal(FARGS, "Can't run multi-threaded MPI simulation yet!");
2481         }
2482         /* and now we just set this; a bit of an ugly hack*/
2483         nnodes = nthreads;
2484     }
2485     /* Check for PME:PP ratio and whether tpr triggers additional output files */
2486     guessPMEratio = inspect_tpr(NFILE, fnm, &rcoulomb);
2487
2488     /* Automatically set -beo options if -eo is set etc. */
2489     couple_files_options(NFILE, fnm);
2490
2491     /* Construct the command line arguments for benchmark runs
2492      * as well as for the simulation run */
2493     if (bThreads)
2494     {
2495         sprintf(bbuf, " -ntmpi %d ", nthreads);
2496     }
2497     else
2498     {
2499         /* This string will be used for MPI runs and will appear after the
2500          * mpirun command. */
2501         if (std::strcmp(procstring[0], "none") != 0)
2502         {
2503             sprintf(bbuf, " -%s %d ", procstring[0], nnodes);
2504         }
2505         else
2506         {
2507             sprintf(bbuf, " ");
2508         }
2509     }
2510
2511     cmd_np = bbuf;
2512
2513     create_command_line_snippets(bAppendFiles,
2514                                  bKeepAndNumCPT,
2515                                  bResetCountersHalfWay,
2516                                  presteps,
2517                                  NFILE,
2518                                  fnm,
2519                                  &cmd_args_bench,
2520                                  &cmd_args_launch,
2521                                  ExtraArgs,
2522                                  deffnm);
2523
2524     /* Prepare to use checkpoint file if requested */
2525     sim_part = 1;
2526     if (opt2bSet("-cpi", NFILE, fnm))
2527     {
2528         const char* filename = opt2fn("-cpi", NFILE, fnm);
2529         int         cpt_sim_part;
2530         read_checkpoint_part_and_step(filename, &cpt_sim_part, &cpt_steps);
2531         if (cpt_sim_part == 0)
2532         {
2533             gmx_fatal(FARGS, "Checkpoint file %s could not be read!", filename);
2534         }
2535         /* tune_pme will run the next part of the simulation */
2536         sim_part = cpt_sim_part + 1;
2537     }
2538
2539     /* Open performance output file and write header info */
2540     fp = gmx_ffopen(opt2fn("-p", NFILE, fnm), "w");
2541
2542     /* Make a quick consistency check of command line parameters */
2543     check_input(nnodes,
2544                 repeats,
2545                 &ntprs,
2546                 &rmin,
2547                 rcoulomb,
2548                 &rmax,
2549                 maxPMEfraction,
2550                 minPMEfraction,
2551                 npme_fixed,
2552                 bench_nsteps,
2553                 fnm,
2554                 NFILE,
2555                 sim_part,
2556                 presteps,
2557                 asize(pa),
2558                 pa);
2559
2560     /* Determine the maximum and minimum number of PME nodes to test,
2561      * the actual list of settings is build in do_the_tests(). */
2562     if ((nnodes > 2) && (npme_fixed < -1))
2563     {
2564         if (0 == std::strcmp(npmevalues_opt[0], "auto"))
2565         {
2566             /* Determine the npme range automatically based on the PME:PP load guess */
2567             if (guessPMEratio > 1.0)
2568             {
2569                 /* More PME than PP work, probably we do not need separate PME nodes at all! */
2570                 maxPMEnodes = nnodes / 2;
2571                 minPMEnodes = nnodes / 2;
2572             }
2573             else
2574             {
2575                 /* PME : PP load is in the range 0..1, let's test around the guess */
2576                 guessPMEnodes = static_cast<int>(nnodes / (1.0 + 1.0 / guessPMEratio));
2577                 minPMEnodes   = static_cast<int>(std::floor(0.7 * guessPMEnodes));
2578                 maxPMEnodes   = static_cast<int>(std::ceil(1.6 * guessPMEnodes));
2579                 maxPMEnodes   = std::min(maxPMEnodes, nnodes / 2);
2580             }
2581         }
2582         else
2583         {
2584             /* Determine the npme range based on user input */
2585             maxPMEnodes = static_cast<int>(std::floor(maxPMEfraction * nnodes));
2586             minPMEnodes = std::max(static_cast<int>(std::floor(minPMEfraction * nnodes)), 0);
2587             fprintf(stdout, "Will try runs with %d ", minPMEnodes);
2588             if (maxPMEnodes != minPMEnodes)
2589             {
2590                 fprintf(stdout, "- %d ", maxPMEnodes);
2591             }
2592             fprintf(stdout,
2593                     "PME-only ranks.\n  Note that the automatic number of PME-only ranks and no "
2594                     "separate PME ranks are always tested.\n");
2595         }
2596     }
2597     else
2598     {
2599         maxPMEnodes = 0;
2600         minPMEnodes = 0;
2601     }
2602
2603     /* Get the commands we need to set up the runs from environment variables */
2604     get_program_paths(bThreads, &cmd_mpirun, &cmd_mdrun);
2605     if (bBenchmark && repeats > 0)
2606     {
2607         check_mdrun_works(bThreads, cmd_mpirun, cmd_np, cmd_mdrun, nullptr != eligible_gpu_ids);
2608     }
2609
2610     /* Print some header info to file */
2611     sep_line(fp);
2612     fprintf(fp, "\n      P E R F O R M A N C E   R E S U L T S\n");
2613     sep_line(fp);
2614     fprintf(fp, "%s for GROMACS %s\n", output_env_get_program_display_name(oenv), gmx_version());
2615     if (!bThreads)
2616     {
2617         fprintf(fp, "Number of ranks         : %d\n", nnodes);
2618         fprintf(fp, "The mpirun command is   : %s\n", cmd_mpirun);
2619         if (std::strcmp(procstring[0], "none") != 0)
2620         {
2621             fprintf(fp, "Passing # of ranks via  : -%s\n", procstring[0]);
2622         }
2623         else
2624         {
2625             fprintf(fp, "Not setting number of ranks in system call\n");
2626         }
2627     }
2628     else
2629     {
2630         fprintf(fp, "Number of threads       : %d\n", nnodes);
2631     }
2632
2633     fprintf(fp, "The mdrun  command is   : %s\n", cmd_mdrun);
2634     fprintf(fp, "mdrun args benchmarks   : %s\n", cmd_args_bench);
2635     fprintf(fp, "Benchmark steps         : ");
2636     fprintf(fp, "%" PRId64, bench_nsteps);
2637     fprintf(fp, "\n");
2638     fprintf(fp, "dlb equilibration steps : %d\n", presteps);
2639     if (sim_part > 1)
2640     {
2641         fprintf(fp, "Checkpoint time step    : ");
2642         fprintf(fp, "%" PRId64, cpt_steps);
2643         fprintf(fp, "\n");
2644     }
2645     fprintf(fp, "mdrun args at launchtime: %s\n", cmd_args_launch);
2646
2647     if (new_sim_nsteps >= 0)
2648     {
2649         bOverwrite = TRUE;
2650         fprintf(stderr, "Note: Simulation input file %s will have ", opt2fn("-so", NFILE, fnm));
2651         fprintf(stderr, "%" PRId64, new_sim_nsteps + cpt_steps);
2652         fprintf(stderr, " steps.\n");
2653         fprintf(fp, "Simulation steps        : ");
2654         fprintf(fp, "%" PRId64, new_sim_nsteps);
2655         fprintf(fp, "\n");
2656     }
2657     if (repeats > 1)
2658     {
2659         fprintf(fp, "Repeats for each test   : %d\n", repeats);
2660     }
2661
2662     if (npme_fixed >= -1)
2663     {
2664         fprintf(fp, "Fixing -npme at         : %d\n", npme_fixed);
2665     }
2666
2667     fprintf(fp, "Input file              : %s\n", opt2fn("-s", NFILE, fnm));
2668     fprintf(fp, "   PME/PP load estimate : %g\n", guessPMEratio);
2669
2670     /* Allocate memory for the inputinfo struct: */
2671     snew(info, 1);
2672     info->nr_inputfiles = ntprs;
2673     for (i = 0; i < ntprs; i++)
2674     {
2675         snew(info->rcoulomb, ntprs);
2676         snew(info->rvdw, ntprs);
2677         snew(info->rlist, ntprs);
2678         snew(info->nkx, ntprs);
2679         snew(info->nky, ntprs);
2680         snew(info->nkz, ntprs);
2681         snew(info->fsx, ntprs);
2682         snew(info->fsy, ntprs);
2683         snew(info->fsz, ntprs);
2684     }
2685     /* Make alternative tpr files to test: */
2686     snew(tpr_names, ntprs);
2687     for (i = 0; i < ntprs; i++)
2688     {
2689         snew(tpr_names[i], STRLEN);
2690     }
2691
2692     /* It can be that ntprs is reduced by make_benchmark_tprs if not enough
2693      * different grids could be found. */
2694     make_benchmark_tprs(opt2fn("-s", NFILE, fnm),
2695                         tpr_names,
2696                         bench_nsteps + presteps,
2697                         cpt_steps,
2698                         rmin,
2699                         rmax,
2700                         bScaleRvdw,
2701                         &ntprs,
2702                         info,
2703                         fp);
2704
2705     /********************************************************************************/
2706     /* Main loop over all scenarios we need to test: tpr files, PME nodes, repeats  */
2707     /********************************************************************************/
2708     snew(perfdata, ntprs);
2709     if (bBenchmark)
2710     {
2711         GMX_RELEASE_ASSERT(npmevalues_opt[0] != nullptr,
2712                            "Options inconsistency; npmevalues_opt[0] is NULL");
2713         do_the_tests(fp,
2714                      tpr_names,
2715                      maxPMEnodes,
2716                      minPMEnodes,
2717                      npme_fixed,
2718                      npmevalues_opt[0],
2719                      perfdata,
2720                      &pmeentries,
2721                      repeats,
2722                      nnodes,
2723                      ntprs,
2724                      bThreads,
2725                      cmd_mpirun,
2726                      cmd_np,
2727                      cmd_mdrun,
2728                      cmd_args_bench,
2729                      fnm,
2730                      NFILE,
2731                      presteps,
2732                      cpt_steps,
2733                      bCheck,
2734                      eligible_gpu_ids);
2735
2736         fprintf(fp, "\nTuning took%8.1f minutes.\n", (gmx_gettime() - seconds) / 60.0);
2737
2738         /* Analyse the results and give a suggestion for optimal settings: */
2739         bKeepTPR = analyze_data(
2740                 fp, opt2fn("-p", NFILE, fnm), perfdata, nnodes, ntprs, pmeentries, repeats, info, &best_tpr, &best_npme);
2741
2742         /* Take the best-performing tpr file and enlarge nsteps to original value */
2743         if (bKeepTPR && !bOverwrite)
2744         {
2745             simulation_tpr = opt2fn("-s", NFILE, fnm);
2746         }
2747         else
2748         {
2749             simulation_tpr = opt2fn("-so", NFILE, fnm);
2750             modify_PMEsettings(bOverwrite ? (new_sim_nsteps + cpt_steps) : info->orig_sim_steps,
2751                                info->orig_init_step,
2752                                tpr_names[best_tpr],
2753                                simulation_tpr);
2754         }
2755
2756         /* Let's get rid of the temporary benchmark input files */
2757         for (i = 0; i < ntprs; i++)
2758         {
2759             fprintf(stdout, "Deleting temporary benchmark input file %s\n", tpr_names[i]);
2760             remove(tpr_names[i]);
2761         }
2762
2763         /* Now start the real simulation if the user requested it ... */
2764         launch_simulation(
2765                 bLaunch, fp, bThreads, cmd_mpirun, cmd_np, cmd_mdrun, cmd_args_launch, simulation_tpr, best_npme, eligible_gpu_ids);
2766     }
2767     gmx_ffclose(fp);
2768
2769     /* ... or simply print the performance results to screen: */
2770     if (!bLaunch)
2771     {
2772         finalize(opt2fn("-p", NFILE, fnm));
2773     }
2774
2775     return 0;
2776 }