src/tools/gmx_tune_pme.c

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