src/gromacs/gmxana/gmx_tune_pme.c

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