src/gromacs/gmxana/gmx_tune_pme.c

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