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 need 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                 if (ecutsVERLET == ir->cutoff_scheme)
1012                 {
1013                     /* With Verlet, the van der Waals radius must always equal the Coulomb radius */
1014                     ir->rvdw = ir->rcoulomb;
1015                 }
1016                 else
1017                 {
1018                     /* For vdw cutoff, rvdw >= rlist */
1019                     ir->rvdw = max(info->rvdw[0], ir->rlist);
1020                 }
1021             }
1022
1023             ir->rlistlong = max_cutoff(ir->rlist, max_cutoff(ir->rvdw, ir->rcoulomb));
1024
1025         } /* end of "if (j != 0)" */
1026
1027         /* for j==0: Save the original settings
1028          * for j >0: Save modified radii and Fourier grids */
1029         info->rcoulomb[j]  = ir->rcoulomb;
1030         info->rvdw[j]      = ir->rvdw;
1031         info->nkx[j]       = ir->nkx;
1032         info->nky[j]       = ir->nky;
1033         info->nkz[j]       = ir->nkz;
1034         info->rlist[j]     = ir->rlist;
1035         info->rlistlong[j] = ir->rlistlong;
1036         info->fsx[j]       = fac*fourierspacing;
1037         info->fsy[j]       = fac*fourierspacing;
1038         info->fsz[j]       = fac*fourierspacing;
1039
1040         /* Write the benchmark tpr file */
1041         strncpy(fn_bench_tprs[j], fn_sim_tpr, strlen(fn_sim_tpr)-strlen(".tpr"));
1042         sprintf(buf, "_bench%.2d.tpr", j);
1043         strcat(fn_bench_tprs[j], buf);
1044         fprintf(stdout, "Writing benchmark tpr %s with nsteps=", fn_bench_tprs[j]);
1045         fprintf(stdout, "%"GMX_PRId64, ir->nsteps);
1046         if (j > 0)
1047         {
1048             fprintf(stdout, ", scaling factor %f\n", fac);
1049         }
1050         else
1051         {
1052             fprintf(stdout, ", unmodified settings\n");
1053         }
1054
1055         write_tpx_state(fn_bench_tprs[j], ir, &state, &mtop);
1056
1057         /* Write information about modified tpr settings to log file */
1058         fprintf(fp, "%4d%10f%10f", j, fac, ir->rcoulomb);
1059         fprintf(fp, "%5d%5d%5d", ir->nkx, ir->nky, ir->nkz);
1060         fprintf(fp, " %9f ", info->fsx[j]);
1061         if (evdwCUT == ir->vdwtype)
1062         {
1063             fprintf(fp, "%10f", ir->rvdw);
1064         }
1065         if (EPME_SWITCHED(ir->coulombtype))
1066         {
1067             fprintf(fp, "%10f", ir->rlist);
1068         }
1069         if (info->rlistlong[0] != max_cutoff(info->rlist[0], max_cutoff(info->rvdw[0], info->rcoulomb[0])) )
1070         {
1071             fprintf(fp, "%10f", ir->rlistlong);
1072         }
1073         fprintf(fp, "  %-14s\n", fn_bench_tprs[j]);
1074
1075         /* Make it clear to the user that some additional settings were modified */
1076         if (!is_equal(ir->rvdw, info->rvdw[0])
1077             || !is_equal(ir->rlistlong, info->rlistlong[0]) )
1078         {
1079             bNote = TRUE;
1080         }
1081     }
1082     if (bNote)
1083     {
1084         fprintf(fp, "\nNote that in addition to the Coulomb radius and the Fourier grid\n"
1085                 "other input settings were also changed (see table above).\n"
1086                 "Please check if the modified settings are appropriate.\n");
1087     }
1088     fflush(stdout);
1089     fflush(fp);
1090     sfree(ir);
1091 }
1092
1093
1094 /* Rename the files we want to keep to some meaningful filename and
1095  * delete the rest */
1096 static void cleanup(const t_filenm *fnm, int nfile, int k, int nnodes,
1097                     int nPMEnodes, int nr, gmx_bool bKeepStderr)
1098 {
1099     char        numstring[STRLEN];
1100     char        newfilename[STRLEN];
1101     const char *fn = NULL;
1102     int         i;
1103     const char *opt;
1104
1105
1106     fprintf(stdout, "Cleaning up, deleting benchmark temp files ...\n");
1107
1108     for (i = 0; i < nfile; i++)
1109     {
1110         opt = (char *)fnm[i].opt;
1111         if (strcmp(opt, "-p") == 0)
1112         {
1113             /* do nothing; keep this file */
1114             ;
1115         }
1116         else if (strcmp(opt, "-bg") == 0)
1117         {
1118             /* Give the log file a nice name so one can later see which parameters were used */
1119             numstring[0] = '\0';
1120             if (nr > 0)
1121             {
1122                 sprintf(numstring, "_%d", nr);
1123             }
1124             sprintf(newfilename, "%s_no%d_np%d_npme%d%s", opt2fn("-bg", nfile, fnm), k, nnodes, nPMEnodes, numstring);
1125             if (gmx_fexist(opt2fn("-bg", nfile, fnm)))
1126             {
1127                 fprintf(stdout, "renaming log file to %s\n", newfilename);
1128                 make_backup(newfilename);
1129                 rename(opt2fn("-bg", nfile, fnm), newfilename);
1130             }
1131         }
1132         else if (strcmp(opt, "-err") == 0)
1133         {
1134             /* This file contains the output of stderr. We want to keep it in
1135              * cases where there have been problems. */
1136             fn           = opt2fn(opt, nfile, fnm);
1137             numstring[0] = '\0';
1138             if (nr > 0)
1139             {
1140                 sprintf(numstring, "_%d", nr);
1141             }
1142             sprintf(newfilename, "%s_no%d_np%d_npme%d%s", fn, k, nnodes, nPMEnodes, numstring);
1143             if (gmx_fexist(fn))
1144             {
1145                 if (bKeepStderr)
1146                 {
1147                     fprintf(stdout, "Saving stderr output in %s\n", newfilename);
1148                     make_backup(newfilename);
1149                     rename(fn, newfilename);
1150                 }
1151                 else
1152                 {
1153                     fprintf(stdout, "Deleting %s\n", fn);
1154                     remove(fn);
1155                 }
1156             }
1157         }
1158         /* Delete the files which are created for each benchmark run: (options -b*) */
1159         else if ( (0 == strncmp(opt, "-b", 2)) && (opt2bSet(opt, nfile, fnm) || !is_optional(&fnm[i])) )
1160         {
1161             fn = opt2fn(opt, nfile, fnm);
1162             if (gmx_fexist(fn))
1163             {
1164                 fprintf(stdout, "Deleting %s\n", fn);
1165                 remove(fn);
1166             }
1167         }
1168     }
1169 }
1170
1171
1172 /* Returns the largest common factor of n1 and n2 */
1173 static int largest_common_factor(int n1, int n2)
1174 {
1175     int factor, nmax;
1176
1177     nmax = min(n1, n2);
1178     for (factor = nmax; factor > 0; factor--)
1179     {
1180         if (0 == (n1 % factor) && 0 == (n2 % factor) )
1181         {
1182             return(factor);
1183         }
1184     }
1185     return 0; /* one for the compiler */
1186 }
1187
1188 enum {
1189     eNpmeAuto, eNpmeAll, eNpmeReduced, eNpmeSubset, eNpmeNr
1190 };
1191
1192 /* Create a list of numbers of PME nodes to test */
1193 static void make_npme_list(
1194         const char *npmevalues_opt, /* Make a complete list with all
1195                                      * possibilities or a short list that keeps only
1196                                      * reasonable numbers of PME nodes                  */
1197         int        *nentries,       /* Number of entries we put in the nPMEnodes list   */
1198         int        *nPMEnodes[],    /* Each entry contains the value for -npme          */
1199         int         nnodes,         /* Total number of nodes to do the tests on         */
1200         int         minPMEnodes,    /* Minimum number of PME nodes                      */
1201         int         maxPMEnodes)    /* Maximum number of PME nodes                      */
1202 {
1203     int i, npme, npp;
1204     int min_factor = 1;   /* We request that npp and npme have this minimal
1205                            * largest common factor (depends on npp)           */
1206     int nlistmax;         /* Max. list size                                   */
1207     int nlist;            /* Actual number of entries in list                 */
1208     int eNPME = 0;
1209
1210
1211     /* Do we need to check all possible values for -npme or is a reduced list enough? */
1212     if (0 == strcmp(npmevalues_opt, "all") )
1213     {
1214         eNPME = eNpmeAll;
1215     }
1216     else if (0 == strcmp(npmevalues_opt, "subset") )
1217     {
1218         eNPME = eNpmeSubset;
1219     }
1220     else /* "auto" or "range" */
1221     {
1222         if (nnodes <= 64)
1223         {
1224             eNPME = eNpmeAll;
1225         }
1226         else if (nnodes < 128)
1227         {
1228             eNPME = eNpmeReduced;
1229         }
1230         else
1231         {
1232             eNPME = eNpmeSubset;
1233         }
1234     }
1235
1236     /* Calculate how many entries we could possibly have (in case of -npme all) */
1237     if (nnodes > 2)
1238     {
1239         nlistmax = maxPMEnodes - minPMEnodes + 3;
1240         if (0 == minPMEnodes)
1241         {
1242             nlistmax--;
1243         }
1244     }
1245     else
1246     {
1247         nlistmax = 1;
1248     }
1249
1250     /* Now make the actual list which is at most of size nlist */
1251     snew(*nPMEnodes, nlistmax);
1252     nlist = 0; /* start counting again, now the real entries in the list */
1253     for (i = 0; i < nlistmax - 2; i++)
1254     {
1255         npme = maxPMEnodes - i;
1256         npp  = nnodes-npme;
1257         switch (eNPME)
1258         {
1259             case eNpmeAll:
1260                 min_factor = 1;
1261                 break;
1262             case eNpmeReduced:
1263                 min_factor = 2;
1264                 break;
1265             case eNpmeSubset:
1266                 /* For 2d PME we want a common largest factor of at least the cube
1267                  * root of the number of PP nodes */
1268                 min_factor = (int) pow(npp, 1.0/3.0);
1269                 break;
1270             default:
1271                 gmx_fatal(FARGS, "Unknown option for eNPME in make_npme_list");
1272                 break;
1273         }
1274         if (largest_common_factor(npp, npme) >= min_factor)
1275         {
1276             (*nPMEnodes)[nlist] = npme;
1277             nlist++;
1278         }
1279     }
1280     /* We always test 0 PME nodes and the automatic number */
1281     *nentries             = nlist + 2;
1282     (*nPMEnodes)[nlist  ] =  0;
1283     (*nPMEnodes)[nlist+1] = -1;
1284
1285     fprintf(stderr, "Will try the following %d different values for -npme:\n", *nentries);
1286     for (i = 0; i < *nentries-1; i++)
1287     {
1288         fprintf(stderr, "%d, ", (*nPMEnodes)[i]);
1289     }
1290     fprintf(stderr, "and %d (auto).\n", (*nPMEnodes)[*nentries-1]);
1291 }
1292
1293
1294 /* Allocate memory to store the performance data */
1295 static void init_perfdata(t_perf *perfdata[], int ntprs, int datasets, int repeats)
1296 {
1297     int i, j, k;
1298
1299
1300     for (k = 0; k < ntprs; k++)
1301     {
1302         snew(perfdata[k], datasets);
1303         for (i = 0; i < datasets; i++)
1304         {
1305             for (j = 0; j < repeats; j++)
1306             {
1307                 snew(perfdata[k][i].Gcycles, repeats);
1308                 snew(perfdata[k][i].ns_per_day, repeats);
1309                 snew(perfdata[k][i].PME_f_load, repeats);
1310             }
1311         }
1312     }
1313 }
1314
1315
1316 /* Check for errors on mdrun -h */
1317 static void make_sure_it_runs(char *mdrun_cmd_line, int length, FILE *fp)
1318 {
1319     char *command, *msg;
1320     int   ret;
1321
1322
1323     snew(command, length +  15);
1324     snew(msg, length + 500);
1325
1326     fprintf(stdout, "Making sure the benchmarks can be executed ...\n");
1327     /* FIXME: mdrun -h no longer actually does anything useful.
1328      * It unconditionally prints the help, ignoring all other options. */
1329     sprintf(command, "%s-h -quiet", mdrun_cmd_line);
1330     ret = gmx_system_call(command);
1331
1332     if (0 != ret)
1333     {
1334         /* To prevent confusion, do not again issue a gmx_fatal here since we already
1335          * get the error message from mdrun itself */
1336         sprintf(msg,    "Cannot run the benchmark simulations! Please check the error message of\n"
1337                 "mdrun for the source of the problem. Did you provide a command line\n"
1338                 "argument that neither g_tune_pme nor mdrun understands? Offending command:\n"
1339                 "\n%s\n\n", command);
1340
1341         fprintf(stderr, "%s", msg);
1342         sep_line(fp);
1343         fprintf(fp, "%s", msg);
1344
1345         exit(ret);
1346     }
1347
1348     sfree(command);
1349     sfree(msg    );
1350 }
1351
1352
1353 static void do_the_tests(
1354         FILE           *fp,             /* General g_tune_pme output file         */
1355         char          **tpr_names,      /* Filenames of the input files to test   */
1356         int             maxPMEnodes,    /* Max fraction of nodes to use for PME   */
1357         int             minPMEnodes,    /* Min fraction of nodes to use for PME   */
1358         int             npme_fixed,     /* If >= -1, test fixed number of PME
1359                                          * nodes only                             */
1360         const char     *npmevalues_opt, /* Which -npme values should be tested    */
1361         t_perf        **perfdata,       /* Here the performace data is stored     */
1362         int            *pmeentries,     /* Entries in the nPMEnodes list          */
1363         int             repeats,        /* Repeat each test this often            */
1364         int             nnodes,         /* Total number of nodes = nPP + nPME     */
1365         int             nr_tprs,        /* Total number of tpr files to test      */
1366         gmx_bool        bThreads,       /* Threads or MPI?                        */
1367         char           *cmd_mpirun,     /* mpirun command string                  */
1368         char           *cmd_np,         /* "-np", "-n", whatever mpirun needs     */
1369         char           *cmd_mdrun,      /* mdrun command string                   */
1370         char           *cmd_args_bench, /* arguments for mdrun in a string        */
1371         const t_filenm *fnm,            /* List of filenames from command line    */
1372         int             nfile,          /* Number of files specified on the cmdl. */
1373         int             presteps,       /* DLB equilibration steps, is checked    */
1374         gmx_int64_t     cpt_steps)      /* Time step counter in the checkpoint    */
1375 {
1376     int      i, nr, k, ret, count = 0, totaltests;
1377     int     *nPMEnodes = NULL;
1378     t_perf  *pd        = NULL;
1379     int      cmdline_length;
1380     char    *command, *cmd_stub;
1381     char     buf[STRLEN];
1382     gmx_bool bResetProblem = FALSE;
1383     gmx_bool bFirst        = TRUE;
1384
1385
1386     /* This string array corresponds to the eParselog enum type at the start
1387      * of this file */
1388     const char* ParseLog[] = {
1389         "OK.",
1390         "Logfile not found!",
1391         "No timings, logfile truncated?",
1392         "Run was terminated.",
1393         "Counters were not reset properly.",
1394         "No DD grid found for these settings.",
1395         "TPX version conflict!",
1396         "mdrun was not started in parallel!",
1397         "Number of PP nodes has a prime factor that is too large.",
1398         "An error occured."
1399     };
1400     char        str_PME_f_load[13];
1401
1402
1403     /* Allocate space for the mdrun command line. 100 extra characters should
1404        be more than enough for the -npme etcetera arguments */
1405     cmdline_length =  strlen(cmd_mpirun)
1406         + strlen(cmd_np)
1407         + strlen(cmd_mdrun)
1408         + strlen(cmd_args_bench)
1409         + strlen(tpr_names[0]) + 100;
1410     snew(command, cmdline_length);
1411     snew(cmd_stub, cmdline_length);
1412
1413     /* Construct the part of the command line that stays the same for all tests: */
1414     if (bThreads)
1415     {
1416         sprintf(cmd_stub, "%s%s", cmd_mdrun, cmd_np);
1417     }
1418     else
1419     {
1420         sprintf(cmd_stub, "%s%s%s ", cmd_mpirun, cmd_np, cmd_mdrun);
1421     }
1422
1423     /* Create a list of numbers of PME nodes to test */
1424     if (npme_fixed < -1)
1425     {
1426         make_npme_list(npmevalues_opt, pmeentries, &nPMEnodes,
1427                        nnodes, minPMEnodes, maxPMEnodes);
1428     }
1429     else
1430     {
1431         *pmeentries  = 1;
1432         snew(nPMEnodes, 1);
1433         nPMEnodes[0] = npme_fixed;
1434         fprintf(stderr, "Will use a fixed number of %d PME-only nodes.\n", nPMEnodes[0]);
1435     }
1436
1437     if (0 == repeats)
1438     {
1439         fprintf(fp, "\nNo benchmarks done since number of repeats (-r) is 0.\n");
1440         gmx_ffclose(fp);
1441         finalize(opt2fn("-p", nfile, fnm));
1442         exit(0);
1443     }
1444
1445     /* Allocate one dataset for each tpr input file: */
1446     init_perfdata(perfdata, nr_tprs, *pmeentries, repeats);
1447
1448     /*****************************************/
1449     /* Main loop over all tpr files to test: */
1450     /*****************************************/
1451     totaltests = nr_tprs*(*pmeentries)*repeats;
1452     for (k = 0; k < nr_tprs; k++)
1453     {
1454         fprintf(fp, "\nIndividual timings for input file %d (%s):\n", k, tpr_names[k]);
1455         fprintf(fp, "PME nodes      Gcycles       ns/day        PME/f    Remark\n");
1456         /* Loop over various numbers of PME nodes: */
1457         for (i = 0; i < *pmeentries; i++)
1458         {
1459             pd = &perfdata[k][i];
1460
1461             /* Loop over the repeats for each scenario: */
1462             for (nr = 0; nr < repeats; nr++)
1463             {
1464                 pd->nPMEnodes = nPMEnodes[i];
1465
1466                 /* Add -npme and -s to the command line and save it. Note that
1467                  * the -passall (if set) options requires cmd_args_bench to be
1468                  * at the end of the command line string */
1469                 snew(pd->mdrun_cmd_line, cmdline_length);
1470                 sprintf(pd->mdrun_cmd_line, "%s-npme %d -s %s %s",
1471                         cmd_stub, pd->nPMEnodes, tpr_names[k], cmd_args_bench);
1472
1473                 /* To prevent that all benchmarks fail due to a show-stopper argument
1474                  * on the mdrun command line, we make a quick check with mdrun -h first */
1475                 if (bFirst)
1476                 {
1477                     make_sure_it_runs(pd->mdrun_cmd_line, cmdline_length, fp);
1478                 }
1479                 bFirst = FALSE;
1480
1481                 /* Do a benchmark simulation: */
1482                 if (repeats > 1)
1483                 {
1484                     sprintf(buf, ", pass %d/%d", nr+1, repeats);
1485                 }
1486                 else
1487                 {
1488                     buf[0] = '\0';
1489                 }
1490                 fprintf(stdout, "\n=== Progress %2.0f%%, tpr %d/%d, run %d/%d%s:\n",
1491                         (100.0*count)/totaltests,
1492                         k+1, nr_tprs, i+1, *pmeentries, buf);
1493                 make_backup(opt2fn("-err", nfile, fnm));
1494                 sprintf(command, "%s 1> /dev/null 2>%s", pd->mdrun_cmd_line, opt2fn("-err", nfile, fnm));
1495                 fprintf(stdout, "%s\n", pd->mdrun_cmd_line);
1496                 gmx_system_call(command);
1497
1498                 /* Collect the performance data from the log file; also check stderr
1499                  * for fatal errors */
1500                 ret = parse_logfile(opt2fn("-bg", nfile, fnm), opt2fn("-err", nfile, fnm),
1501                                     pd, nr, presteps, cpt_steps, nnodes);
1502                 if ((presteps > 0) && (ret == eParselogResetProblem))
1503                 {
1504                     bResetProblem = TRUE;
1505                 }
1506
1507                 if (-1 == pd->nPMEnodes)
1508                 {
1509                     sprintf(buf, "(%3d)", pd->guessPME);
1510                 }
1511                 else
1512                 {
1513                     sprintf(buf, "     ");
1514                 }
1515
1516                 /* Nicer output */
1517                 if (pd->PME_f_load[nr] > 0.0)
1518                 {
1519                     sprintf(str_PME_f_load, "%12.3f", pd->PME_f_load[nr]);
1520                 }
1521                 else
1522                 {
1523                     sprintf(str_PME_f_load, "%s", "         -  ");
1524                 }
1525
1526                 /* Write the data we got to disk */
1527                 fprintf(fp, "%4d%s %12.3f %12.3f %s    %s", pd->nPMEnodes,
1528                         buf, pd->Gcycles[nr], pd->ns_per_day[nr], str_PME_f_load, ParseLog[ret]);
1529                 if (!(ret == eParselogOK || ret == eParselogNoDDGrid || ret == eParselogNotFound) )
1530                 {
1531                     fprintf(fp, " Check %s file for problems.", ret == eParselogFatal ? "err" : "log");
1532                 }
1533                 fprintf(fp, "\n");
1534                 fflush(fp);
1535                 count++;
1536
1537                 /* Do some cleaning up and delete the files we do not need any more */
1538                 cleanup(fnm, nfile, k, nnodes, pd->nPMEnodes, nr, ret == eParselogFatal);
1539
1540                 /* If the first run with this number of processors already failed, do not try again: */
1541                 if (pd->Gcycles[0] <= 0.0 && repeats > 1)
1542                 {
1543                     fprintf(stdout, "Skipping remaining passes of unsuccessful setting, see log file for details.\n");
1544                     count += repeats-(nr+1);
1545                     break;
1546                 }
1547             } /* end of repeats loop */
1548         }     /* end of -npme loop */
1549     }         /* end of tpr file loop */
1550
1551     if (bResetProblem)
1552     {
1553         sep_line(fp);
1554         fprintf(fp, "WARNING: The cycle and time step counters could not be reset properly. ");
1555         sep_line(fp);
1556     }
1557     sfree(command);
1558     sfree(cmd_stub);
1559 }
1560
1561
1562 static void check_input(
1563         int             nnodes,
1564         int             repeats,
1565         int            *ntprs,
1566         real           *rmin,
1567         real            rcoulomb,
1568         real           *rmax,
1569         real            maxPMEfraction,
1570         real            minPMEfraction,
1571         int             npme_fixed,
1572         gmx_int64_t     bench_nsteps,
1573         const t_filenm *fnm,
1574         int             nfile,
1575         int             sim_part,
1576         int             presteps,
1577         int             npargs,
1578         t_pargs        *pa)
1579 {
1580     int old;
1581
1582
1583     /* Make sure the input file exists */
1584     if (!gmx_fexist(opt2fn("-s", nfile, fnm)))
1585     {
1586         gmx_fatal(FARGS, "File %s not found.", opt2fn("-s", nfile, fnm));
1587     }
1588
1589     /* Make sure that the checkpoint file is not overwritten during benchmarking */
1590     if ( (0 == strcmp(opt2fn("-cpi", nfile, fnm), opt2fn("-bcpo", nfile, fnm)) ) && (sim_part > 1) )
1591     {
1592         gmx_fatal(FARGS, "Checkpoint input (-cpi) and benchmark checkpoint output (-bcpo) files must not be identical.\n"
1593                   "The checkpoint input file must not be overwritten during the benchmarks.\n");
1594     }
1595
1596     /* Make sure that repeats is >= 0 (if == 0, only write tpr files) */
1597     if (repeats < 0)
1598     {
1599         gmx_fatal(FARGS, "Number of repeats < 0!");
1600     }
1601
1602     /* Check number of nodes */
1603     if (nnodes < 1)
1604     {
1605         gmx_fatal(FARGS, "Number of nodes/threads must be a positive integer.");
1606     }
1607
1608     /* Automatically choose -ntpr if not set */
1609     if (*ntprs < 1)
1610     {
1611         if (nnodes < 16)
1612         {
1613             *ntprs = 1;
1614         }
1615         else
1616         {
1617             *ntprs = 3;
1618             /* Set a reasonable scaling factor for rcoulomb */
1619             if (*rmax <= 0)
1620             {
1621                 *rmax = rcoulomb * 1.2;
1622             }
1623         }
1624         fprintf(stderr, "Will test %d tpr file%s.\n", *ntprs, *ntprs == 1 ? "" : "s");
1625     }
1626     else
1627     {
1628         if (1 == *ntprs)
1629         {
1630             fprintf(stderr, "Note: Choose ntpr>1 to shift PME load between real and reciprocal space.\n");
1631         }
1632     }
1633
1634     /* Make shure that rmin <= rcoulomb <= rmax */
1635     if (*rmin <= 0)
1636     {
1637         *rmin = rcoulomb;
1638     }
1639     if (*rmax <= 0)
1640     {
1641         *rmax = rcoulomb;
1642     }
1643     if (!(*rmin <= *rmax) )
1644     {
1645         gmx_fatal(FARGS, "Please choose the Coulomb radii such that rmin <= rmax.\n"
1646                   "rmin = %g, rmax = %g, actual rcoul from .tpr file = %g\n", *rmin, *rmax, rcoulomb);
1647     }
1648     /* Add test scenarios if rmin or rmax were set */
1649     if (*ntprs <= 2)
1650     {
1651         if (!is_equal(*rmin, rcoulomb) && (*ntprs == 1) )
1652         {
1653             (*ntprs)++;
1654             fprintf(stderr, "NOTE: Setting -rmin to %g changed -ntpr to %d\n",
1655                     *rmin, *ntprs);
1656         }
1657         if (!is_equal(*rmax, rcoulomb) && (*ntprs == 1) )
1658         {
1659             (*ntprs)++;
1660             fprintf(stderr, "NOTE: Setting -rmax to %g changed -ntpr to %d\n",
1661                     *rmax, *ntprs);
1662         }
1663     }
1664     old = *ntprs;
1665     /* If one of rmin, rmax is set, we need 2 tpr files at minimum */
1666     if (!is_equal(*rmax, rcoulomb) || !is_equal(*rmin, rcoulomb) )
1667     {
1668         *ntprs = max(*ntprs, 2);
1669     }
1670
1671     /* If both rmin, rmax are set, we need 3 tpr files at minimum */
1672     if (!is_equal(*rmax, rcoulomb) && !is_equal(*rmin, rcoulomb) )
1673     {
1674         *ntprs = max(*ntprs, 3);
1675     }
1676
1677     if (old != *ntprs)
1678     {
1679         fprintf(stderr, "NOTE: Your rmin, rmax setting changed -ntpr to %d\n", *ntprs);
1680     }
1681
1682     if (*ntprs > 1)
1683     {
1684         if (is_equal(*rmin, rcoulomb) && is_equal(rcoulomb, *rmax)) /* We have just a single rc */
1685         {
1686             fprintf(stderr, "WARNING: Resetting -ntpr to 1 since no Coulomb radius scaling is requested.\n"
1687                     "Please set rmin < rmax to test Coulomb radii in the [rmin, rmax] interval\n"
1688                     "with correspondingly adjusted PME grid settings\n");
1689             *ntprs = 1;
1690         }
1691     }
1692
1693     /* Check whether max and min fraction are within required values */
1694     if (maxPMEfraction > 0.5 || maxPMEfraction < 0)
1695     {
1696         gmx_fatal(FARGS, "-max must be between 0 and 0.5");
1697     }
1698     if (minPMEfraction > 0.5 || minPMEfraction < 0)
1699     {
1700         gmx_fatal(FARGS, "-min must be between 0 and 0.5");
1701     }
1702     if (maxPMEfraction < minPMEfraction)
1703     {
1704         gmx_fatal(FARGS, "-max must be larger or equal to -min");
1705     }
1706
1707     /* Check whether the number of steps - if it was set - has a reasonable value */
1708     if (bench_nsteps < 0)
1709     {
1710         gmx_fatal(FARGS, "Number of steps must be positive.");
1711     }
1712
1713     if (bench_nsteps > 10000 || bench_nsteps < 100)
1714     {
1715         fprintf(stderr, "WARNING: steps=");
1716         fprintf(stderr, "%"GMX_PRId64, bench_nsteps);
1717         fprintf(stderr, ". Are you sure you want to perform so %s steps for each benchmark?\n", (bench_nsteps < 100) ? "few" : "many");
1718     }
1719
1720     if (presteps < 0)
1721     {
1722         gmx_fatal(FARGS, "Cannot have a negative number of presteps.\n");
1723     }
1724
1725     /* Check for rcoulomb scaling if more than one .tpr file is tested */
1726     if (*ntprs > 1)
1727     {
1728         if (*rmin/rcoulomb < 0.75 || *rmax/rcoulomb > 1.25)
1729         {
1730             fprintf(stderr, "WARNING: Applying extreme scaling factor. I hope you know what you are doing.\n");
1731         }
1732     }
1733
1734     /* If a fixed number of PME nodes is set we do rcoulomb and PME gird tuning
1735      * only. We need to check whether the requested number of PME-only nodes
1736      * makes sense. */
1737     if (npme_fixed > -1)
1738     {
1739         /* No more than 50% of all nodes can be assigned as PME-only nodes. */
1740         if (2*npme_fixed > nnodes)
1741         {
1742             gmx_fatal(FARGS, "Cannot have more than %d PME-only nodes for a total of %d nodes (you chose %d).\n",
1743                       nnodes/2, nnodes, npme_fixed);
1744         }
1745         if ((npme_fixed > 0) && (5*npme_fixed < nnodes))
1746         {
1747             fprintf(stderr, "WARNING: Only %g percent of the nodes are assigned as PME-only nodes.\n",
1748                     100.0*((real)npme_fixed / (real)nnodes));
1749         }
1750         if (opt2parg_bSet("-min", npargs, pa) || opt2parg_bSet("-max", npargs, pa))
1751         {
1752             fprintf(stderr, "NOTE: The -min, -max, and -npme options have no effect when a\n"
1753                     "      fixed number of PME-only nodes is requested with -fix.\n");
1754         }
1755     }
1756 }
1757
1758
1759 /* Returns TRUE when "opt" is needed at launch time */
1760 static gmx_bool is_launch_file(char *opt, gmx_bool bSet)
1761 {
1762     if (0 == strncmp(opt, "-swap", 5))
1763     {
1764         return bSet;
1765     }
1766
1767     /* Apart from the input .tpr and the output log files we need all options that
1768      * were set on the command line and that do not start with -b */
1769     if    (0 == strncmp(opt, "-b", 2) || 0 == strncmp(opt, "-s", 2)
1770            || 0 == strncmp(opt, "-err", 4) || 0 == strncmp(opt, "-p", 2) )
1771     {
1772         return FALSE;
1773     }
1774
1775     return bSet;
1776 }
1777
1778
1779 /* Returns TRUE when "opt" defines a file which is needed for the benchmarks runs */
1780 static gmx_bool is_bench_file(char *opt, gmx_bool bSet, gmx_bool bOptional, gmx_bool bIsOutput)
1781 {
1782     /* Apart from the input .tpr, all files starting with "-b" are for
1783      * _b_enchmark files exclusively */
1784     if (0 == strncmp(opt, "-s", 2))
1785     {
1786         return FALSE;
1787     }
1788
1789     if (0 == strncmp(opt, "-b", 2) || 0 == strncmp(opt, "-s", 2))
1790     {
1791         if (!bOptional || bSet)
1792         {
1793             return TRUE;
1794         }
1795         else
1796         {
1797             return FALSE;
1798         }
1799     }
1800     else
1801     {
1802         if (bIsOutput)
1803         {
1804             return FALSE;
1805         }
1806         else
1807         {
1808             if (bSet) /* These are additional input files like -cpi -ei */
1809             {
1810                 return TRUE;
1811             }
1812             else
1813             {
1814                 return FALSE;
1815             }
1816         }
1817     }
1818 }
1819
1820
1821 /* Adds 'buf' to 'str' */
1822 static void add_to_string(char **str, char *buf)
1823 {
1824     int len;
1825
1826
1827     len = strlen(*str) + strlen(buf) + 1;
1828     srenew(*str, len);
1829     strcat(*str, buf);
1830 }
1831
1832
1833 /* Create the command line for the benchmark as well as for the real run */
1834 static void create_command_line_snippets(
1835         gmx_bool  bAppendFiles,
1836         gmx_bool  bKeepAndNumCPT,
1837         gmx_bool  bResetHWay,
1838         int       presteps,
1839         int       nfile,
1840         t_filenm  fnm[],
1841         char     *cmd_args_bench[],  /* command line arguments for benchmark runs */
1842         char     *cmd_args_launch[], /* command line arguments for simulation run */
1843         char      extra_args[])      /* Add this to the end of the command line */
1844 {
1845     int         i;
1846     char       *opt;
1847     const char *name;
1848     char        strbuf[STRLEN];
1849
1850
1851     /* strlen needs at least '\0' as a string: */
1852     snew(*cmd_args_bench, 1);
1853     snew(*cmd_args_launch, 1);
1854     *cmd_args_launch[0] = '\0';
1855     *cmd_args_bench[0]  = '\0';
1856
1857
1858     /*******************************************/
1859     /* 1. Process other command line arguments */
1860     /*******************************************/
1861     if (presteps > 0)
1862     {
1863         /* Add equilibration steps to benchmark options */
1864         sprintf(strbuf, "-resetstep %d ", presteps);
1865         add_to_string(cmd_args_bench, strbuf);
1866     }
1867     /* These switches take effect only at launch time */
1868     if (FALSE == bAppendFiles)
1869     {
1870         add_to_string(cmd_args_launch, "-noappend ");
1871     }
1872     if (bKeepAndNumCPT)
1873     {
1874         add_to_string(cmd_args_launch, "-cpnum ");
1875     }
1876     if (bResetHWay)
1877     {
1878         add_to_string(cmd_args_launch, "-resethway ");
1879     }
1880
1881     /********************/
1882     /* 2. Process files */
1883     /********************/
1884     for (i = 0; i < nfile; i++)
1885     {
1886         opt  = (char *)fnm[i].opt;
1887         name = opt2fn(opt, nfile, fnm);
1888
1889         /* Strbuf contains the options, now let's sort out where we need that */
1890         sprintf(strbuf, "%s %s ", opt, name);
1891
1892         if (is_bench_file(opt, opt2bSet(opt, nfile, fnm), is_optional(&fnm[i]), is_output(&fnm[i])) )
1893         {
1894             /* All options starting with -b* need the 'b' removed,
1895              * therefore overwrite strbuf */
1896             if (0 == strncmp(opt, "-b", 2))
1897             {
1898                 sprintf(strbuf, "-%s %s ", &opt[2], name);
1899             }
1900
1901             add_to_string(cmd_args_bench, strbuf);
1902         }
1903
1904         if (is_launch_file(opt, opt2bSet(opt, nfile, fnm)) )
1905         {
1906             add_to_string(cmd_args_launch, strbuf);
1907         }
1908     }
1909
1910     add_to_string(cmd_args_bench, extra_args);
1911     add_to_string(cmd_args_launch, extra_args);
1912 }
1913
1914
1915 /* Set option opt */
1916 static void setopt(const char *opt, int nfile, t_filenm fnm[])
1917 {
1918     int i;
1919
1920     for (i = 0; (i < nfile); i++)
1921     {
1922         if (strcmp(opt, fnm[i].opt) == 0)
1923         {
1924             fnm[i].flag |= ffSET;
1925         }
1926     }
1927 }
1928
1929
1930 /* This routine inspects the tpr file and ...
1931  * 1. checks for output files that get triggered by a tpr option. These output
1932  *    files are marked as 'set' to allow for a proper cleanup after each
1933  *    tuning run.
1934  * 2. returns the PME:PP load ratio
1935  * 3. returns rcoulomb from the tpr */
1936 static float inspect_tpr(int nfile, t_filenm fnm[], real *rcoulomb)
1937 {
1938     gmx_bool     bPull;     /* Is pulling requested in .tpr file?             */
1939     gmx_bool     bTpi;      /* Is test particle insertion requested?          */
1940     gmx_bool     bFree;     /* Is a free energy simulation requested?         */
1941     gmx_bool     bNM;       /* Is a normal mode analysis requested?           */
1942     gmx_bool     bSwap;     /* Is water/ion position swapping requested?      */
1943     t_inputrec   ir;
1944     t_state      state;
1945     gmx_mtop_t   mtop;
1946
1947
1948     /* Check tpr file for options that trigger extra output files */
1949     read_tpx_state(opt2fn("-s", nfile, fnm), &ir, &state, NULL, &mtop);
1950     bPull = (epullNO != ir.ePull);
1951     bFree = (efepNO  != ir.efep );
1952     bNM   = (eiNM    == ir.eI   );
1953     bSwap = (eswapNO != ir.eSwapCoords);
1954     bTpi  = EI_TPI(ir.eI);
1955
1956     /* Set these output files on the tuning command-line */
1957     if (bPull)
1958     {
1959         setopt("-pf", nfile, fnm);
1960         setopt("-px", nfile, fnm);
1961     }
1962     if (bFree)
1963     {
1964         setopt("-dhdl", nfile, fnm);
1965     }
1966     if (bTpi)
1967     {
1968         setopt("-tpi", nfile, fnm);
1969         setopt("-tpid", nfile, fnm);
1970     }
1971     if (bNM)
1972     {
1973         setopt("-mtx", nfile, fnm);
1974     }
1975     if (bSwap)
1976     {
1977         setopt("-swap", nfile, fnm);
1978     }
1979
1980     *rcoulomb = ir.rcoulomb;
1981
1982     /* Return the estimate for the number of PME nodes */
1983     return pme_load_estimate(&mtop, &ir, state.box);
1984 }
1985
1986
1987 static void couple_files_options(int nfile, t_filenm fnm[])
1988 {
1989     int      i;
1990     gmx_bool bSet, bBench;
1991     char    *opt;
1992     char     buf[20];
1993
1994
1995     for (i = 0; i < nfile; i++)
1996     {
1997         opt    = (char *)fnm[i].opt;
1998         bSet   = ((fnm[i].flag & ffSET) != 0);
1999         bBench = (0 == strncmp(opt, "-b", 2));
2000
2001         /* Check optional files */
2002         /* If e.g. -eo is set, then -beo also needs to be set */
2003         if (is_optional(&fnm[i]) && bSet && !bBench)
2004         {
2005             sprintf(buf, "-b%s", &opt[1]);
2006             setopt(buf, nfile, fnm);
2007         }
2008         /* If -beo is set, then -eo also needs to be! */
2009         if (is_optional(&fnm[i]) && bSet && bBench)
2010         {
2011             sprintf(buf, "-%s", &opt[2]);
2012             setopt(buf, nfile, fnm);
2013         }
2014     }
2015 }
2016
2017
2018 #define BENCHSTEPS (1000)
2019
2020 int gmx_tune_pme(int argc, char *argv[])
2021 {
2022     const char     *desc[] = {
2023         "For a given number [TT]-np[tt] or [TT]-ntmpi[tt] of processors/threads, [THISMODULE] systematically",
2024         "times [gmx-mdrun] with various numbers of PME-only nodes and determines",
2025         "which setting is fastest. It will also test whether performance can",
2026         "be enhanced by shifting load from the reciprocal to the real space",
2027         "part of the Ewald sum. ",
2028         "Simply pass your [TT].tpr[tt] file to [THISMODULE] together with other options",
2029         "for [gmx-mdrun] as needed.[PAR]",
2030         "Which executables are used can be set in the environment variables",
2031         "MPIRUN and MDRUN. If these are not present, 'mpirun' and 'mdrun'",
2032         "will be used as defaults. Note that for certain MPI frameworks you",
2033         "need to provide a machine- or hostfile. This can also be passed",
2034         "via the MPIRUN variable, e.g.[PAR]",
2035         "[TT]export MPIRUN=\"/usr/local/mpirun -machinefile hosts\"[tt][PAR]",
2036         "Please call [THISMODULE] with the normal options you would pass to",
2037         "[gmx-mdrun] and add [TT]-np[tt] for the number of processors to perform the",
2038         "tests on, or [TT]-ntmpi[tt] for the number of threads. You can also add [TT]-r[tt]",
2039         "to repeat each test several times to get better statistics. [PAR]",
2040         "[THISMODULE] can test various real space / reciprocal space workloads",
2041         "for you. With [TT]-ntpr[tt] you control how many extra [TT].tpr[tt] files will be",
2042         "written with enlarged cutoffs and smaller Fourier grids respectively.",
2043         "Typically, the first test (number 0) will be with the settings from the input",
2044         "[TT].tpr[tt] file; the last test (number [TT]ntpr[tt]) will have the Coulomb cutoff",
2045         "specified by [TT]-rmax[tt] with a somwhat smaller PME grid at the same time. ",
2046         "In this last test, the Fourier spacing is multiplied with [TT]rmax[tt]/rcoulomb. ",
2047         "The remaining [TT].tpr[tt] files will have equally-spaced Coulomb radii (and Fourier "
2048         "spacings) between these extremes. [BB]Note[bb] that you can set [TT]-ntpr[tt] to 1",
2049         "if you just seek the optimal number of PME-only nodes; in that case",
2050         "your input [TT].tpr[tt] file will remain unchanged.[PAR]",
2051         "For the benchmark runs, the default of 1000 time steps should suffice for most",
2052         "MD systems. The dynamic load balancing needs about 100 time steps",
2053         "to adapt to local load imbalances, therefore the time step counters",
2054         "are by default reset after 100 steps. For large systems (>1M atoms), as well as ",
2055         "for a higher accuarcy of the measurements, you should set [TT]-resetstep[tt] to a higher value.",
2056         "From the 'DD' load imbalance entries in the md.log output file you",
2057         "can tell after how many steps the load is sufficiently balanced. Example call:[PAR]"
2058         "[TT]gmx tune_pme -np 64 -s protein.tpr -launch[tt][PAR]",
2059         "After calling [gmx-mdrun] several times, detailed performance information",
2060         "is available in the output file [TT]perf.out[tt].",
2061         "[BB]Note[bb] that during the benchmarks, a couple of temporary files are written",
2062         "(options [TT]-b[tt]*), these will be automatically deleted after each test.[PAR]",
2063         "If you want the simulation to be started automatically with the",
2064         "optimized parameters, use the command line option [TT]-launch[tt].[PAR]",
2065     };
2066
2067     int             nnodes         = 1;
2068     int             repeats        = 2;
2069     int             pmeentries     = 0; /* How many values for -npme do we actually test for each tpr file */
2070     real            maxPMEfraction = 0.50;
2071     real            minPMEfraction = 0.25;
2072     int             maxPMEnodes, minPMEnodes;
2073     float           guessPMEratio;                    /* guessed PME:PP ratio based on the tpr file */
2074     float           guessPMEnodes;
2075     int             npme_fixed     = -2;              /* If >= -1, use only this number
2076                                                        * of PME-only nodes                */
2077     int             ntprs          = 0;
2078     real            rmin           = 0.0, rmax = 0.0; /* min and max value for rcoulomb if scaling is requested */
2079     real            rcoulomb       = -1.0;            /* Coulomb radius as set in .tpr file */
2080     gmx_bool        bScaleRvdw     = TRUE;
2081     gmx_int64_t     bench_nsteps   = BENCHSTEPS;
2082     gmx_int64_t     new_sim_nsteps = -1;  /* -1 indicates: not set by the user */
2083     gmx_int64_t     cpt_steps      = 0;   /* Step counter in .cpt input file   */
2084     int             presteps       = 100; /* Do a full cycle reset after presteps steps */
2085     gmx_bool        bOverwrite     = FALSE, bKeepTPR;
2086     gmx_bool        bLaunch        = FALSE;
2087     char           *ExtraArgs      = NULL;
2088     char          **tpr_names      = NULL;
2089     const char     *simulation_tpr = NULL;
2090     int             best_npme, best_tpr;
2091     int             sim_part = 1; /* For benchmarks with checkpoint files */
2092     char            bbuf[STRLEN];
2093
2094     /* Default program names if nothing else is found */
2095     char         *cmd_mpirun = NULL, *cmd_mdrun = NULL;
2096     char         *cmd_args_bench, *cmd_args_launch;
2097     char         *cmd_np = NULL;
2098
2099     t_perf      **perfdata = NULL;
2100     t_inputinfo  *info;
2101     int           i;
2102     FILE         *fp;
2103     t_commrec    *cr;
2104
2105     /* Print out how long the tuning took */
2106     double          seconds;
2107
2108     static t_filenm fnm[] = {
2109         /* g_tune_pme */
2110         { efOUT, "-p",      "perf",     ffWRITE },
2111         { efLOG, "-err",    "bencherr", ffWRITE },
2112         { efTPX, "-so",     "tuned",    ffWRITE },
2113         /* mdrun: */
2114         { efTPX, NULL,      NULL,       ffREAD },
2115         { efTRN, "-o",      NULL,       ffWRITE },
2116         { efCOMPRESSED, "-x", NULL,     ffOPTWR },
2117         { efCPT, "-cpi",    NULL,       ffOPTRD },
2118         { efCPT, "-cpo",    NULL,       ffOPTWR },
2119         { efSTO, "-c",      "confout",  ffWRITE },
2120         { efEDR, "-e",      "ener",     ffWRITE },
2121         { efLOG, "-g",      "md",       ffWRITE },
2122         { efXVG, "-dhdl",   "dhdl",     ffOPTWR },
2123         { efXVG, "-field",  "field",    ffOPTWR },
2124         { efXVG, "-table",  "table",    ffOPTRD },
2125         { efXVG, "-tabletf", "tabletf",   ffOPTRD },
2126         { efXVG, "-tablep", "tablep",   ffOPTRD },
2127         { efXVG, "-tableb", "table",    ffOPTRD },
2128         { efTRX, "-rerun",  "rerun",    ffOPTRD },
2129         { efXVG, "-tpi",    "tpi",      ffOPTWR },
2130         { efXVG, "-tpid",   "tpidist",  ffOPTWR },
2131         { efEDI, "-ei",     "sam",      ffOPTRD },
2132         { efXVG, "-eo",     "edsam",    ffOPTWR },
2133         { efXVG, "-devout", "deviatie", ffOPTWR },
2134         { efXVG, "-runav",  "runaver",  ffOPTWR },
2135         { efXVG, "-px",     "pullx",    ffOPTWR },
2136         { efXVG, "-pf",     "pullf",    ffOPTWR },
2137         { efXVG, "-ro",     "rotation", ffOPTWR },
2138         { efLOG, "-ra",     "rotangles", ffOPTWR },
2139         { efLOG, "-rs",     "rotslabs", ffOPTWR },
2140         { efLOG, "-rt",     "rottorque", ffOPTWR },
2141         { efMTX, "-mtx",    "nm",       ffOPTWR },
2142         { efNDX, "-dn",     "dipole",   ffOPTWR },
2143         { efXVG, "-swap",   "swapions", ffOPTWR },
2144         /* Output files that are deleted after each benchmark run */
2145         { efTRN, "-bo",     "bench",    ffWRITE },
2146         { efXTC, "-bx",     "bench",    ffWRITE },
2147         { efCPT, "-bcpo",   "bench",    ffWRITE },
2148         { efSTO, "-bc",     "bench",    ffWRITE },
2149         { efEDR, "-be",     "bench",    ffWRITE },
2150         { efLOG, "-bg",     "bench",    ffWRITE },
2151         { efXVG, "-beo",    "benchedo", ffOPTWR },
2152         { efXVG, "-bdhdl",  "benchdhdl", ffOPTWR },
2153         { efXVG, "-bfield", "benchfld", ffOPTWR },
2154         { efXVG, "-btpi",   "benchtpi", ffOPTWR },
2155         { efXVG, "-btpid",  "benchtpid", ffOPTWR },
2156         { efXVG, "-bdevout", "benchdev", ffOPTWR },
2157         { efXVG, "-brunav", "benchrnav", ffOPTWR },
2158         { efXVG, "-bpx",    "benchpx",  ffOPTWR },
2159         { efXVG, "-bpf",    "benchpf",  ffOPTWR },
2160         { efXVG, "-bro",    "benchrot", ffOPTWR },
2161         { efLOG, "-bra",    "benchrota", ffOPTWR },
2162         { efLOG, "-brs",    "benchrots", ffOPTWR },
2163         { efLOG, "-brt",    "benchrott", ffOPTWR },
2164         { efMTX, "-bmtx",   "benchn",   ffOPTWR },
2165         { efNDX, "-bdn",    "bench",    ffOPTWR },
2166         { efXVG, "-bswap",  "benchswp", ffOPTWR }
2167     };
2168
2169     gmx_bool        bThreads     = FALSE;
2170
2171     int             nthreads = 1;
2172
2173     const char     *procstring[] =
2174     { NULL, "-np", "-n", "none", NULL };
2175     const char     *npmevalues_opt[] =
2176     { NULL, "auto", "all", "subset", NULL };
2177
2178     gmx_bool     bAppendFiles          = TRUE;
2179     gmx_bool     bKeepAndNumCPT        = FALSE;
2180     gmx_bool     bResetCountersHalfWay = FALSE;
2181     gmx_bool     bBenchmark            = TRUE;
2182
2183     output_env_t oenv = NULL;
2184
2185     t_pargs      pa[] = {
2186         /***********************/
2187         /* g_tune_pme options: */
2188         /***********************/
2189         { "-np",       FALSE, etINT,  {&nnodes},
2190           "Number of nodes to run the tests on (must be > 2 for separate PME nodes)" },
2191         { "-npstring", FALSE, etENUM, {procstring},
2192           "Specify the number of processors to [TT]$MPIRUN[tt] using this string"},
2193         { "-ntmpi",    FALSE, etINT,  {&nthreads},
2194           "Number of MPI-threads to run the tests on (turns MPI & mpirun off)"},
2195         { "-r",        FALSE, etINT,  {&repeats},
2196           "Repeat each test this often" },
2197         { "-max",      FALSE, etREAL, {&maxPMEfraction},
2198           "Max fraction of PME nodes to test with" },
2199         { "-min",      FALSE, etREAL, {&minPMEfraction},
2200           "Min fraction of PME nodes to test with" },
2201         { "-npme",     FALSE, etENUM, {npmevalues_opt},
2202           "Within -min and -max, benchmark all possible values for [TT]-npme[tt], or just a reasonable subset. "
2203           "Auto neglects -min and -max and chooses reasonable values around a guess for npme derived from the .tpr"},
2204         { "-fix",      FALSE, etINT,  {&npme_fixed},
2205           "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."},
2206         { "-rmax",     FALSE, etREAL, {&rmax},
2207           "If >0, maximal rcoulomb for -ntpr>1 (rcoulomb upscaling results in fourier grid downscaling)" },
2208         { "-rmin",     FALSE, etREAL, {&rmin},
2209           "If >0, minimal rcoulomb for -ntpr>1" },
2210         { "-scalevdw",  FALSE, etBOOL, {&bScaleRvdw},
2211           "Scale rvdw along with rcoulomb"},
2212         { "-ntpr",     FALSE, etINT,  {&ntprs},
2213           "Number of [TT].tpr[tt] files to benchmark. Create this many files with different rcoulomb scaling factors depending on -rmin and -rmax. "
2214           "If < 1, automatically choose the number of [TT].tpr[tt] files to test" },
2215         { "-steps",    FALSE, etINT64, {&bench_nsteps},
2216           "Take timings for this many steps in the benchmark runs" },
2217         { "-resetstep", FALSE, etINT,  {&presteps},
2218           "Let dlb equilibrate this many steps before timings are taken (reset cycle counters after this many steps)" },
2219         { "-simsteps", FALSE, etINT64, {&new_sim_nsteps},
2220           "If non-negative, perform this many steps in the real run (overwrites nsteps from [TT].tpr[tt], add [TT].cpt[tt] steps)" },
2221         { "-launch",   FALSE, etBOOL, {&bLaunch},
2222           "Launch the real simulation after optimization" },
2223         { "-bench",    FALSE, etBOOL, {&bBenchmark},
2224           "Run the benchmarks or just create the input [TT].tpr[tt] files?" },
2225         /******************/
2226         /* mdrun options: */
2227         /******************/
2228         /* We let g_tune_pme parse and understand these options, because we need to
2229          * prevent that they appear on the mdrun command line for the benchmarks */
2230         { "-append",   FALSE, etBOOL, {&bAppendFiles},
2231           "Append to previous output files when continuing from checkpoint instead of adding the simulation part number to all file names (for launch only)" },
2232         { "-cpnum",    FALSE, etBOOL, {&bKeepAndNumCPT},
2233           "Keep and number checkpoint files (launch only)" },
2234         { "-resethway", FALSE, etBOOL, {&bResetCountersHalfWay},
2235           "HIDDENReset the cycle counters after half the number of steps or halfway [TT]-maxh[tt] (launch only)" }
2236     };
2237
2238 #define NFILE asize(fnm)
2239
2240     seconds = gmx_gettime();
2241
2242     if (!parse_common_args(&argc, argv, PCA_NOEXIT_ON_ARGS,
2243                            NFILE, fnm, asize(pa), pa, asize(desc), desc,
2244                            0, NULL, &oenv))
2245     {
2246         return 0;
2247     }
2248
2249     /* Store the remaining unparsed command line entries in a string which
2250      * is then attached to the mdrun command line */
2251     snew(ExtraArgs, 1);
2252     ExtraArgs[0] = '\0';
2253     for (i = 1; i < argc; i++) /* argc will now be 1 if everything was understood */
2254     {
2255         add_to_string(&ExtraArgs, argv[i]);
2256         add_to_string(&ExtraArgs, " ");
2257     }
2258
2259     if (opt2parg_bSet("-ntmpi", asize(pa), pa))
2260     {
2261         bThreads = TRUE;
2262         if (opt2parg_bSet("-npstring", asize(pa), pa))
2263         {
2264             fprintf(stderr, "WARNING: -npstring has no effect when using threads.\n");
2265         }
2266
2267         if (nnodes > 1)
2268         {
2269             gmx_fatal(FARGS, "Can't run multi-threaded MPI simulation yet!");
2270         }
2271         /* and now we just set this; a bit of an ugly hack*/
2272         nnodes = nthreads;
2273     }
2274     /* Check for PME:PP ratio and whether tpr triggers additional output files */
2275     guessPMEratio = inspect_tpr(NFILE, fnm, &rcoulomb);
2276
2277     /* Automatically set -beo options if -eo is set etc. */
2278     couple_files_options(NFILE, fnm);
2279
2280     /* Construct the command line arguments for benchmark runs
2281      * as well as for the simulation run */
2282     if (bThreads)
2283     {
2284         sprintf(bbuf, " -ntmpi %d ", nthreads);
2285     }
2286     else
2287     {
2288         /* This string will be used for MPI runs and will appear after the
2289          * mpirun command. */
2290         if (strcmp(procstring[0], "none") != 0)
2291         {
2292             sprintf(bbuf, " %s %d ", procstring[0], nnodes);
2293         }
2294         else
2295         {
2296             sprintf(bbuf, " ");
2297         }
2298     }
2299
2300     cmd_np = bbuf;
2301
2302     create_command_line_snippets(bAppendFiles, bKeepAndNumCPT, bResetCountersHalfWay, presteps,
2303                                  NFILE, fnm, &cmd_args_bench, &cmd_args_launch, ExtraArgs);
2304
2305     /* Read in checkpoint file if requested */
2306     sim_part = 1;
2307     if (opt2bSet("-cpi", NFILE, fnm))
2308     {
2309         snew(cr, 1);
2310         cr->duty = DUTY_PP; /* makes the following routine happy */
2311         read_checkpoint_simulation_part(opt2fn("-cpi", NFILE, fnm),
2312                                         &sim_part, &cpt_steps, cr,
2313                                         FALSE, NFILE, fnm, NULL, NULL);
2314         sfree(cr);
2315         sim_part++;
2316         /* sim_part will now be 1 if no checkpoint file was found */
2317         if (sim_part <= 1)
2318         {
2319             gmx_fatal(FARGS, "Checkpoint file %s not found!", opt2fn("-cpi", NFILE, fnm));
2320         }
2321     }
2322
2323     /* Open performance output file and write header info */
2324     fp = gmx_ffopen(opt2fn("-p", NFILE, fnm), "w");
2325
2326     /* Make a quick consistency check of command line parameters */
2327     check_input(nnodes, repeats, &ntprs, &rmin, rcoulomb, &rmax,
2328                 maxPMEfraction, minPMEfraction, npme_fixed,
2329                 bench_nsteps, fnm, NFILE, sim_part, presteps,
2330                 asize(pa), pa);
2331
2332     /* Determine the maximum and minimum number of PME nodes to test,
2333      * the actual list of settings is build in do_the_tests(). */
2334     if ((nnodes > 2) && (npme_fixed < -1))
2335     {
2336         if (0 == strcmp(npmevalues_opt[0], "auto"))
2337         {
2338             /* Determine the npme range automatically based on the PME:PP load guess */
2339             if (guessPMEratio > 1.0)
2340             {
2341                 /* More PME than PP work, probably we do not need separate PME nodes at all! */
2342                 maxPMEnodes = nnodes/2;
2343                 minPMEnodes = nnodes/2;
2344             }
2345             else
2346             {
2347                 /* PME : PP load is in the range 0..1, let's test around the guess */
2348                 guessPMEnodes = nnodes/(1.0 + 1.0/guessPMEratio);
2349                 minPMEnodes   = floor(0.7*guessPMEnodes);
2350                 maxPMEnodes   =  ceil(1.6*guessPMEnodes);
2351                 maxPMEnodes   = min(maxPMEnodes, nnodes/2);
2352             }
2353         }
2354         else
2355         {
2356             /* Determine the npme range based on user input */
2357             maxPMEnodes = floor(maxPMEfraction*nnodes);
2358             minPMEnodes = max(floor(minPMEfraction*nnodes), 0);
2359             fprintf(stdout, "Will try runs with %d ", minPMEnodes);
2360             if (maxPMEnodes != minPMEnodes)
2361             {
2362                 fprintf(stdout, "- %d ", maxPMEnodes);
2363             }
2364             fprintf(stdout, "PME-only nodes.\n  Note that the automatic number of PME-only nodes and no separate PME nodes are always tested.\n");
2365         }
2366     }
2367     else
2368     {
2369         maxPMEnodes = 0;
2370         minPMEnodes = 0;
2371     }
2372
2373     /* Get the commands we need to set up the runs from environment variables */
2374     get_program_paths(bThreads, &cmd_mpirun, &cmd_mdrun);
2375     if (bBenchmark && repeats > 0)
2376     {
2377         check_mdrun_works(bThreads, cmd_mpirun, cmd_np, cmd_mdrun);
2378     }
2379
2380     /* Print some header info to file */
2381     sep_line(fp);
2382     fprintf(fp, "\n      P E R F O R M A N C E   R E S U L T S\n");
2383     sep_line(fp);
2384     fprintf(fp, "%s for Gromacs %s\n", ShortProgram(), GromacsVersion());
2385     if (!bThreads)
2386     {
2387         fprintf(fp, "Number of nodes         : %d\n", nnodes);
2388         fprintf(fp, "The mpirun command is   : %s\n", cmd_mpirun);
2389         if (strcmp(procstring[0], "none") != 0)
2390         {
2391             fprintf(fp, "Passing # of nodes via  : %s\n", procstring[0]);
2392         }
2393         else
2394         {
2395             fprintf(fp, "Not setting number of nodes in system call\n");
2396         }
2397     }
2398     else
2399     {
2400         fprintf(fp, "Number of threads       : %d\n", nnodes);
2401     }
2402
2403     fprintf(fp, "The mdrun  command is   : %s\n", cmd_mdrun);
2404     fprintf(fp, "mdrun args benchmarks   : %s\n", cmd_args_bench);
2405     fprintf(fp, "Benchmark steps         : ");
2406     fprintf(fp, "%"GMX_PRId64, bench_nsteps);
2407     fprintf(fp, "\n");
2408     fprintf(fp, "dlb equilibration steps : %d\n", presteps);
2409     if (sim_part > 1)
2410     {
2411         fprintf(fp, "Checkpoint time step    : ");
2412         fprintf(fp, "%"GMX_PRId64, cpt_steps);
2413         fprintf(fp, "\n");
2414     }
2415     fprintf(fp, "mdrun args at launchtime: %s\n", cmd_args_launch);
2416
2417     if (new_sim_nsteps >= 0)
2418     {
2419         bOverwrite = TRUE;
2420         fprintf(stderr, "Note: Simulation input file %s will have ", opt2fn("-so", NFILE, fnm));
2421         fprintf(stderr, "%"GMX_PRId64, new_sim_nsteps+cpt_steps);
2422         fprintf(stderr, " steps.\n");
2423         fprintf(fp, "Simulation steps        : ");
2424         fprintf(fp, "%"GMX_PRId64, new_sim_nsteps);
2425         fprintf(fp, "\n");
2426     }
2427     if (repeats > 1)
2428     {
2429         fprintf(fp, "Repeats for each test   : %d\n", repeats);
2430     }
2431
2432     if (npme_fixed >= -1)
2433     {
2434         fprintf(fp, "Fixing -npme at         : %d\n", npme_fixed);
2435     }
2436
2437     fprintf(fp, "Input file              : %s\n", opt2fn("-s", NFILE, fnm));
2438     fprintf(fp, "   PME/PP load estimate : %g\n", guessPMEratio);
2439
2440     /* Allocate memory for the inputinfo struct: */
2441     snew(info, 1);
2442     info->nr_inputfiles = ntprs;
2443     for (i = 0; i < ntprs; i++)
2444     {
2445         snew(info->rcoulomb, ntprs);
2446         snew(info->rvdw, ntprs);
2447         snew(info->rlist, ntprs);
2448         snew(info->rlistlong, ntprs);
2449         snew(info->nkx, ntprs);
2450         snew(info->nky, ntprs);
2451         snew(info->nkz, ntprs);
2452         snew(info->fsx, ntprs);
2453         snew(info->fsy, ntprs);
2454         snew(info->fsz, ntprs);
2455     }
2456     /* Make alternative tpr files to test: */
2457     snew(tpr_names, ntprs);
2458     for (i = 0; i < ntprs; i++)
2459     {
2460         snew(tpr_names[i], STRLEN);
2461     }
2462
2463     /* It can be that ntprs is reduced by make_benchmark_tprs if not enough
2464      * different grids could be found. */
2465     make_benchmark_tprs(opt2fn("-s", NFILE, fnm), tpr_names, bench_nsteps+presteps,
2466                         cpt_steps, rmin, rmax, bScaleRvdw, &ntprs, info, fp);
2467
2468     /********************************************************************************/
2469     /* Main loop over all scenarios we need to test: tpr files, PME nodes, repeats  */
2470     /********************************************************************************/
2471     snew(perfdata, ntprs);
2472     if (bBenchmark)
2473     {
2474         do_the_tests(fp, tpr_names, maxPMEnodes, minPMEnodes, npme_fixed, npmevalues_opt[0], perfdata, &pmeentries,
2475                      repeats, nnodes, ntprs, bThreads, cmd_mpirun, cmd_np, cmd_mdrun,
2476                      cmd_args_bench, fnm, NFILE, presteps, cpt_steps);
2477
2478         fprintf(fp, "\nTuning took%8.1f minutes.\n", (gmx_gettime()-seconds)/60.0);
2479
2480         /* Analyse the results and give a suggestion for optimal settings: */
2481         bKeepTPR = analyze_data(fp, opt2fn("-p", NFILE, fnm), perfdata, nnodes, ntprs, pmeentries,
2482                                 repeats, info, &best_tpr, &best_npme);
2483
2484         /* Take the best-performing tpr file and enlarge nsteps to original value */
2485         if (bKeepTPR && !bOverwrite)
2486         {
2487             simulation_tpr = opt2fn("-s", NFILE, fnm);
2488         }
2489         else
2490         {
2491             simulation_tpr = opt2fn("-so", NFILE, fnm);
2492             modify_PMEsettings(bOverwrite ? (new_sim_nsteps+cpt_steps) : info->orig_sim_steps,
2493                                info->orig_init_step, tpr_names[best_tpr], simulation_tpr);
2494         }
2495
2496         /* Let's get rid of the temporary benchmark input files */
2497         for (i = 0; i < ntprs; i++)
2498         {
2499             fprintf(stdout, "Deleting temporary benchmark input file %s\n", tpr_names[i]);
2500             remove(tpr_names[i]);
2501         }
2502
2503         /* Now start the real simulation if the user requested it ... */
2504         launch_simulation(bLaunch, fp, bThreads, cmd_mpirun, cmd_np, cmd_mdrun,
2505                           cmd_args_launch, simulation_tpr, best_npme);
2506     }
2507     gmx_ffclose(fp);
2508
2509     /* ... or simply print the performance results to screen: */
2510     if (!bLaunch)
2511     {
2512         finalize(opt2fn("-p", NFILE, fnm));
2513     }
2514
2515     return 0;
2516 }