src/tools/gmx_tune_pme.c

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