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