src/tools/gmx_tune_pme.c

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