src/gromacs/gmxana/gmx_tune_pme.c

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