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