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
  42 #ifdef HAVE_SYS_TIME_H
  43 #include <sys/time.h>
  44 #endif
  45
  46 #include "gromacs/commandline/pargs.h"
  47 #include "gromacs/fileio/tpxio.h"
  48 #include "gromacs/gmxana/gmx_ana.h"
  49 #include "gromacs/legacyheaders/calcgrid.h"
  50 #include "gromacs/legacyheaders/checkpoint.h"
  51 #include "gromacs/legacyheaders/inputrec.h"
  52 #include "gromacs/legacyheaders/macros.h"
  53 #include "gromacs/legacyheaders/names.h"
  54 #include "gromacs/legacyheaders/perf_est.h"
  55 #include "gromacs/legacyheaders/readinp.h"
  56 #include "gromacs/legacyheaders/typedefs.h"
  57 #include "gromacs/legacyheaders/types/commrec.h"
  58 #include "gromacs/math/utilities.h"
  59 #include "gromacs/math/vec.h"
  60 #include "gromacs/timing/walltime_accounting.h"
  61 #include "gromacs/utility/baseversion.h"
  62 #include "gromacs/utility/cstringutil.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,
1312                 "Cannot run the first benchmark simulation! Please check the error message of\n"
1313                 "mdrun for the source of the problem. Did you provide a command line\n"
1314                 "argument that neither gmx tune_pme nor mdrun understands? If you're\n"
1315                 "sure your command line should work, you can bypass this check with \n"
1316                 "gmx tune_pme -nocheck. The failing command was:\n"
1317                 "\n%s\n\n", command);
1318
1319         fprintf(stderr, "%s", msg);
1320         sep_line(fp);
1321         fprintf(fp, "%s", msg);
1322
1323         exit(ret);
1324     }
1325     fprintf(stdout, "Benchmarks can be executed!\n");
1326
1327     /* Clean up the benchmark output files we just created */
1328     fprintf(stdout, "Cleaning up ...\n");
1329     remove_if_exists(opt2fn("-bc", nfile, fnm));
1330     remove_if_exists(opt2fn("-be", nfile, fnm));
1331     remove_if_exists(opt2fn("-bcpo", nfile, fnm));
1332     remove_if_exists(opt2fn("-bg", nfile, fnm));
1333
1334     sfree(command);
1335     sfree(msg    );
1336 }
1337
1338
1339 static void do_the_tests(
1340         FILE           *fp,             /* General g_tune_pme output file         */
1341         char          **tpr_names,      /* Filenames of the input files to test   */
1342         int             maxPMEnodes,    /* Max fraction of nodes to use for PME   */
1343         int             minPMEnodes,    /* Min fraction of nodes to use for PME   */
1344         int             npme_fixed,     /* If >= -1, test fixed number of PME
1345                                          * nodes only                             */
1346         const char     *npmevalues_opt, /* Which -npme values should be tested    */
1347         t_perf        **perfdata,       /* Here the performance data is stored    */
1348         int            *pmeentries,     /* Entries in the nPMEnodes list          */
1349         int             repeats,        /* Repeat each test this often            */
1350         int             nnodes,         /* Total number of nodes = nPP + nPME     */
1351         int             nr_tprs,        /* Total number of tpr files to test      */
1352         gmx_bool        bThreads,       /* Threads or MPI?                        */
1353         char           *cmd_mpirun,     /* mpirun command string                  */
1354         char           *cmd_np,         /* "-np", "-n", whatever mpirun needs     */
1355         char           *cmd_mdrun,      /* mdrun command string                   */
1356         char           *cmd_args_bench, /* arguments for mdrun in a string        */
1357         const t_filenm *fnm,            /* List of filenames from command line    */
1358         int             nfile,          /* Number of files specified on the cmdl. */
1359         int             presteps,       /* DLB equilibration steps, is checked    */
1360         gmx_int64_t     cpt_steps,      /* Time step counter in the checkpoint    */
1361         gmx_bool        bCheck)         /* Check whether benchmark mdrun works    */
1362 {
1363     int      i, nr, k, ret, count = 0, totaltests;
1364     int     *nPMEnodes = NULL;
1365     t_perf  *pd        = NULL;
1366     int      cmdline_length;
1367     char    *command, *cmd_stub;
1368     char     buf[STRLEN];
1369     gmx_bool bResetProblem = FALSE;
1370     gmx_bool bFirst        = TRUE;
1371
1372
1373     /* This string array corresponds to the eParselog enum type at the start
1374      * of this file */
1375     const char* ParseLog[] = {
1376         "OK.",
1377         "Logfile not found!",
1378         "No timings, logfile truncated?",
1379         "Run was terminated.",
1380         "Counters were not reset properly.",
1381         "No DD grid found for these settings.",
1382         "TPX version conflict!",
1383         "mdrun was not started in parallel!",
1384         "Number of PP ranks has a prime factor that is too large.",
1385         "An error occured."
1386     };
1387     char        str_PME_f_load[13];
1388
1389
1390     /* Allocate space for the mdrun command line. 100 extra characters should
1391        be more than enough for the -npme etcetera arguments */
1392     cmdline_length =  strlen(cmd_mpirun)
1393         + strlen(cmd_np)
1394         + strlen(cmd_mdrun)
1395         + strlen(cmd_args_bench)
1396         + strlen(tpr_names[0]) + 100;
1397     snew(command, cmdline_length);
1398     snew(cmd_stub, cmdline_length);
1399
1400     /* Construct the part of the command line that stays the same for all tests: */
1401     if (bThreads)
1402     {
1403         sprintf(cmd_stub, "%s%s", cmd_mdrun, cmd_np);
1404     }
1405     else
1406     {
1407         sprintf(cmd_stub, "%s%s%s ", cmd_mpirun, cmd_np, cmd_mdrun);
1408     }
1409
1410     /* Create a list of numbers of PME nodes to test */
1411     if (npme_fixed < -1)
1412     {
1413         make_npme_list(npmevalues_opt, pmeentries, &nPMEnodes,
1414                        nnodes, minPMEnodes, maxPMEnodes);
1415     }
1416     else
1417     {
1418         *pmeentries  = 1;
1419         snew(nPMEnodes, 1);
1420         nPMEnodes[0] = npme_fixed;
1421         fprintf(stderr, "Will use a fixed number of %d PME-only ranks.\n", nPMEnodes[0]);
1422     }
1423
1424     if (0 == repeats)
1425     {
1426         fprintf(fp, "\nNo benchmarks done since number of repeats (-r) is 0.\n");
1427         gmx_ffclose(fp);
1428         finalize(opt2fn("-p", nfile, fnm));
1429         exit(0);
1430     }
1431
1432     /* Allocate one dataset for each tpr input file: */
1433     init_perfdata(perfdata, nr_tprs, *pmeentries, repeats);
1434
1435     /*****************************************/
1436     /* Main loop over all tpr files to test: */
1437     /*****************************************/
1438     totaltests = nr_tprs*(*pmeentries)*repeats;
1439     for (k = 0; k < nr_tprs; k++)
1440     {
1441         fprintf(fp, "\nIndividual timings for input file %d (%s):\n", k, tpr_names[k]);
1442         fprintf(fp, "PME ranks      Gcycles       ns/day        PME/f    Remark\n");
1443         /* Loop over various numbers of PME nodes: */
1444         for (i = 0; i < *pmeentries; i++)
1445         {
1446             pd = &perfdata[k][i];
1447
1448             /* Loop over the repeats for each scenario: */
1449             for (nr = 0; nr < repeats; nr++)
1450             {
1451                 pd->nPMEnodes = nPMEnodes[i];
1452
1453                 /* Add -npme and -s to the command line and save it. Note that
1454                  * the -passall (if set) options requires cmd_args_bench to be
1455                  * at the end of the command line string */
1456                 snew(pd->mdrun_cmd_line, cmdline_length);
1457                 sprintf(pd->mdrun_cmd_line, "%s-npme %d -s %s %s",
1458                         cmd_stub, pd->nPMEnodes, tpr_names[k], cmd_args_bench);
1459
1460                 /* To prevent that all benchmarks fail due to a show-stopper argument
1461                  * on the mdrun command line, we make a quick check first.
1462                  * This check can be turned off in cases where the automatically chosen
1463                  * number of PME-only ranks leads to a number of PP ranks for which no
1464                  * decomposition can be found (e.g. for large prime numbers) */
1465                 if (bFirst && bCheck)
1466                 {
1467                     make_sure_it_runs(pd->mdrun_cmd_line, cmdline_length, fp, fnm, nfile);
1468                 }
1469                 bFirst = FALSE;
1470
1471                 /* Do a benchmark simulation: */
1472                 if (repeats > 1)
1473                 {
1474                     sprintf(buf, ", pass %d/%d", nr+1, repeats);
1475                 }
1476                 else
1477                 {
1478                     buf[0] = '\0';
1479                 }
1480                 fprintf(stdout, "\n=== Progress %2.0f%%, tpr %d/%d, run %d/%d%s:\n",
1481                         (100.0*count)/totaltests,
1482                         k+1, nr_tprs, i+1, *pmeentries, buf);
1483                 make_backup(opt2fn("-err", nfile, fnm));
1484                 sprintf(command, "%s 1> /dev/null 2>%s", pd->mdrun_cmd_line, opt2fn("-err", nfile, fnm));
1485                 fprintf(stdout, "%s\n", pd->mdrun_cmd_line);
1486                 gmx_system_call(command);
1487
1488                 /* Collect the performance data from the log file; also check stderr
1489                  * for fatal errors */
1490                 ret = parse_logfile(opt2fn("-bg", nfile, fnm), opt2fn("-err", nfile, fnm),
1491                                     pd, nr, presteps, cpt_steps, nnodes);
1492                 if ((presteps > 0) && (ret == eParselogResetProblem))
1493                 {
1494                     bResetProblem = TRUE;
1495                 }
1496
1497                 if (-1 == pd->nPMEnodes)
1498                 {
1499                     sprintf(buf, "(%3d)", pd->guessPME);
1500                 }
1501                 else
1502                 {
1503                     sprintf(buf, "     ");
1504                 }
1505
1506                 /* Nicer output */
1507                 if (pd->PME_f_load[nr] > 0.0)
1508                 {
1509                     sprintf(str_PME_f_load, "%12.3f", pd->PME_f_load[nr]);
1510                 }
1511                 else
1512                 {
1513                     sprintf(str_PME_f_load, "%s", "         -  ");
1514                 }
1515
1516                 /* Write the data we got to disk */
1517                 fprintf(fp, "%4d%s %12.3f %12.3f %s    %s", pd->nPMEnodes,
1518                         buf, pd->Gcycles[nr], pd->ns_per_day[nr], str_PME_f_load, ParseLog[ret]);
1519                 if (!(ret == eParselogOK || ret == eParselogNoDDGrid || ret == eParselogNotFound) )
1520                 {
1521                     fprintf(fp, " Check %s file for problems.", ret == eParselogFatal ? "err" : "log");
1522                 }
1523                 fprintf(fp, "\n");
1524                 fflush(fp);
1525                 count++;
1526
1527                 /* Do some cleaning up and delete the files we do not need any more */
1528                 cleanup(fnm, nfile, k, nnodes, pd->nPMEnodes, nr, ret == eParselogFatal);
1529
1530                 /* If the first run with this number of processors already failed, do not try again: */
1531                 if (pd->Gcycles[0] <= 0.0 && repeats > 1)
1532                 {
1533                     fprintf(stdout, "Skipping remaining passes of unsuccessful setting, see log file for details.\n");
1534                     count += repeats-(nr+1);
1535                     break;
1536                 }
1537             } /* end of repeats loop */
1538         }     /* end of -npme loop */
1539     }         /* end of tpr file loop */
1540
1541     if (bResetProblem)
1542     {
1543         sep_line(fp);
1544         fprintf(fp, "WARNING: The cycle and time step counters could not be reset properly. ");
1545         sep_line(fp);
1546     }
1547     sfree(command);
1548     sfree(cmd_stub);
1549 }
1550
1551
1552 static void check_input(
1553         int             nnodes,
1554         int             repeats,
1555         int            *ntprs,
1556         real           *rmin,
1557         real            rcoulomb,
1558         real           *rmax,
1559         real            maxPMEfraction,
1560         real            minPMEfraction,
1561         int             npme_fixed,
1562         gmx_int64_t     bench_nsteps,
1563         const t_filenm *fnm,
1564         int             nfile,
1565         int             sim_part,
1566         int             presteps,
1567         int             npargs,
1568         t_pargs        *pa)
1569 {
1570     int old;
1571
1572
1573     /* Make sure the input file exists */
1574     if (!gmx_fexist(opt2fn("-s", nfile, fnm)))
1575     {
1576         gmx_fatal(FARGS, "File %s not found.", opt2fn("-s", nfile, fnm));
1577     }
1578
1579     /* Make sure that the checkpoint file is not overwritten during benchmarking */
1580     if ( (0 == strcmp(opt2fn("-cpi", nfile, fnm), opt2fn("-bcpo", nfile, fnm)) ) && (sim_part > 1) )
1581     {
1582         gmx_fatal(FARGS, "Checkpoint input (-cpi) and benchmark checkpoint output (-bcpo) files must not be identical.\n"
1583                   "The checkpoint input file must not be overwritten during the benchmarks.\n");
1584     }
1585
1586     /* Make sure that repeats is >= 0 (if == 0, only write tpr files) */
1587     if (repeats < 0)
1588     {
1589         gmx_fatal(FARGS, "Number of repeats < 0!");
1590     }
1591
1592     /* Check number of nodes */
1593     if (nnodes < 1)
1594     {
1595         gmx_fatal(FARGS, "Number of ranks/threads must be a positive integer.");
1596     }
1597
1598     /* Automatically choose -ntpr if not set */
1599     if (*ntprs < 1)
1600     {
1601         if (nnodes < 16)
1602         {
1603             *ntprs = 1;
1604         }
1605         else
1606         {
1607             *ntprs = 3;
1608             /* Set a reasonable scaling factor for rcoulomb */
1609             if (*rmax <= 0)
1610             {
1611                 *rmax = rcoulomb * 1.2;
1612             }
1613         }
1614         fprintf(stderr, "Will test %d tpr file%s.\n", *ntprs, *ntprs == 1 ? "" : "s");
1615     }
1616     else
1617     {
1618         if (1 == *ntprs)
1619         {
1620             fprintf(stderr, "Note: Choose ntpr>1 to shift PME load between real and reciprocal space.\n");
1621         }
1622     }
1623
1624     /* Make shure that rmin <= rcoulomb <= rmax */
1625     if (*rmin <= 0)
1626     {
1627         *rmin = rcoulomb;
1628     }
1629     if (*rmax <= 0)
1630     {
1631         *rmax = rcoulomb;
1632     }
1633     if (!(*rmin <= *rmax) )
1634     {
1635         gmx_fatal(FARGS, "Please choose the Coulomb radii such that rmin <= rmax.\n"
1636                   "rmin = %g, rmax = %g, actual rcoul from .tpr file = %g\n", *rmin, *rmax, rcoulomb);
1637     }
1638     /* Add test scenarios if rmin or rmax were set */
1639     if (*ntprs <= 2)
1640     {
1641         if (!gmx_within_tol(*rmin, rcoulomb, GMX_REAL_EPS) && (*ntprs == 1) )
1642         {
1643             (*ntprs)++;
1644             fprintf(stderr, "NOTE: Setting -rmin to %g changed -ntpr to %d\n",
1645                     *rmin, *ntprs);
1646         }
1647         if (!gmx_within_tol(*rmax, rcoulomb, GMX_REAL_EPS) && (*ntprs == 1) )
1648         {
1649             (*ntprs)++;
1650             fprintf(stderr, "NOTE: Setting -rmax to %g changed -ntpr to %d\n",
1651                     *rmax, *ntprs);
1652         }
1653     }
1654     old = *ntprs;
1655     /* If one of rmin, rmax is set, we need 2 tpr files at minimum */
1656     if (!gmx_within_tol(*rmax, rcoulomb, GMX_REAL_EPS) || !gmx_within_tol(*rmin, rcoulomb, GMX_REAL_EPS) )
1657     {
1658         *ntprs = max(*ntprs, 2);
1659     }
1660
1661     /* If both rmin, rmax are set, we need 3 tpr files at minimum */
1662     if (!gmx_within_tol(*rmax, rcoulomb, GMX_REAL_EPS) && !gmx_within_tol(*rmin, rcoulomb, GMX_REAL_EPS) )
1663     {
1664         *ntprs = max(*ntprs, 3);
1665     }
1666
1667     if (old != *ntprs)
1668     {
1669         fprintf(stderr, "NOTE: Your rmin, rmax setting changed -ntpr to %d\n", *ntprs);
1670     }
1671
1672     if (*ntprs > 1)
1673     {
1674         if (gmx_within_tol(*rmin, rcoulomb, GMX_REAL_EPS) && gmx_within_tol(rcoulomb, *rmax, GMX_REAL_EPS)) /* We have just a single rc */
1675         {
1676             fprintf(stderr, "WARNING: Resetting -ntpr to 1 since no Coulomb radius scaling is requested.\n"
1677                     "Please set rmin < rmax to test Coulomb radii in the [rmin, rmax] interval\n"
1678                     "with correspondingly adjusted PME grid settings\n");
1679             *ntprs = 1;
1680         }
1681     }
1682
1683     /* Check whether max and min fraction are within required values */
1684     if (maxPMEfraction > 0.5 || maxPMEfraction < 0)
1685     {
1686         gmx_fatal(FARGS, "-max must be between 0 and 0.5");
1687     }
1688     if (minPMEfraction > 0.5 || minPMEfraction < 0)
1689     {
1690         gmx_fatal(FARGS, "-min must be between 0 and 0.5");
1691     }
1692     if (maxPMEfraction < minPMEfraction)
1693     {
1694         gmx_fatal(FARGS, "-max must be larger or equal to -min");
1695     }
1696
1697     /* Check whether the number of steps - if it was set - has a reasonable value */
1698     if (bench_nsteps < 0)
1699     {
1700         gmx_fatal(FARGS, "Number of steps must be positive.");
1701     }
1702
1703     if (bench_nsteps > 10000 || bench_nsteps < 100)
1704     {
1705         fprintf(stderr, "WARNING: steps=");
1706         fprintf(stderr, "%"GMX_PRId64, bench_nsteps);
1707         fprintf(stderr, ". Are you sure you want to perform so %s steps for each benchmark?\n", (bench_nsteps < 100) ? "few" : "many");
1708     }
1709
1710     if (presteps < 0)
1711     {
1712         gmx_fatal(FARGS, "Cannot have a negative number of presteps.\n");
1713     }
1714
1715     /* Check for rcoulomb scaling if more than one .tpr file is tested */
1716     if (*ntprs > 1)
1717     {
1718         if (*rmin/rcoulomb < 0.75 || *rmax/rcoulomb > 1.25)
1719         {
1720             fprintf(stderr, "WARNING: Applying extreme scaling factor. I hope you know what you are doing.\n");
1721         }
1722     }
1723
1724     /* If a fixed number of PME nodes is set we do rcoulomb and PME gird tuning
1725      * only. We need to check whether the requested number of PME-only nodes
1726      * makes sense. */
1727     if (npme_fixed > -1)
1728     {
1729         /* No more than 50% of all nodes can be assigned as PME-only nodes. */
1730         if (2*npme_fixed > nnodes)
1731         {
1732             gmx_fatal(FARGS, "Cannot have more than %d PME-only ranks for a total of %d ranks (you chose %d).\n",
1733                       nnodes/2, nnodes, npme_fixed);
1734         }
1735         if ((npme_fixed > 0) && (5*npme_fixed < nnodes))
1736         {
1737             fprintf(stderr, "WARNING: Only %g percent of the ranks are assigned as PME-only ranks.\n",
1738                     100.0*((real)npme_fixed / (real)nnodes));
1739         }
1740         if (opt2parg_bSet("-min", npargs, pa) || opt2parg_bSet("-max", npargs, pa))
1741         {
1742             fprintf(stderr, "NOTE: The -min, -max, and -npme options have no effect when a\n"
1743                     "      fixed number of PME-only ranks is requested with -fix.\n");
1744         }
1745     }
1746 }
1747
1748
1749 /* Returns TRUE when "opt" is needed at launch time */
1750 static gmx_bool is_launch_file(char *opt, gmx_bool bSet)
1751 {
1752     if (0 == strncmp(opt, "-swap", 5))
1753     {
1754         return bSet;
1755     }
1756
1757     /* Apart from the input .tpr and the output log files we need all options that
1758      * were set on the command line and that do not start with -b */
1759     if    (0 == strncmp(opt, "-b", 2) || 0 == strncmp(opt, "-s", 2)
1760            || 0 == strncmp(opt, "-err", 4) || 0 == strncmp(opt, "-p", 2) )
1761     {
1762         return FALSE;
1763     }
1764
1765     return bSet;
1766 }
1767
1768
1769 /* Returns TRUE when "opt" defines a file which is needed for the benchmarks runs */
1770 static gmx_bool is_bench_file(char *opt, gmx_bool bSet, gmx_bool bOptional, gmx_bool bIsOutput)
1771 {
1772     /* Apart from the input .tpr, all files starting with "-b" are for
1773      * _b_enchmark files exclusively */
1774     if (0 == strncmp(opt, "-s", 2))
1775     {
1776         return FALSE;
1777     }
1778
1779     if (0 == strncmp(opt, "-b", 2) || 0 == strncmp(opt, "-s", 2))
1780     {
1781         if (!bOptional || bSet)
1782         {
1783             return TRUE;
1784         }
1785         else
1786         {
1787             return FALSE;
1788         }
1789     }
1790     else
1791     {
1792         if (bIsOutput)
1793         {
1794             return FALSE;
1795         }
1796         else
1797         {
1798             if (bSet) /* These are additional input files like -cpi -ei */
1799             {
1800                 return TRUE;
1801             }
1802             else
1803             {
1804                 return FALSE;
1805             }
1806         }
1807     }
1808 }
1809
1810
1811 /* Adds 'buf' to 'str' */
1812 static void add_to_string(char **str, char *buf)
1813 {
1814     int len;
1815
1816
1817     len = strlen(*str) + strlen(buf) + 1;
1818     srenew(*str, len);
1819     strcat(*str, buf);
1820 }
1821
1822
1823 /* Create the command line for the benchmark as well as for the real run */
1824 static void create_command_line_snippets(
1825         gmx_bool  bAppendFiles,
1826         gmx_bool  bKeepAndNumCPT,
1827         gmx_bool  bResetHWay,
1828         int       presteps,
1829         int       nfile,
1830         t_filenm  fnm[],
1831         char     *cmd_args_bench[],  /* command line arguments for benchmark runs */
1832         char     *cmd_args_launch[], /* command line arguments for simulation run */
1833         char      extra_args[])      /* Add this to the end of the command line */
1834 {
1835     int         i;
1836     char       *opt;
1837     const char *name;
1838     char        strbuf[STRLEN];
1839
1840
1841     /* strlen needs at least '\0' as a string: */
1842     snew(*cmd_args_bench, 1);
1843     snew(*cmd_args_launch, 1);
1844     *cmd_args_launch[0] = '\0';
1845     *cmd_args_bench[0]  = '\0';
1846
1847
1848     /*******************************************/
1849     /* 1. Process other command line arguments */
1850     /*******************************************/
1851     if (presteps > 0)
1852     {
1853         /* Add equilibration steps to benchmark options */
1854         sprintf(strbuf, "-resetstep %d ", presteps);
1855         add_to_string(cmd_args_bench, strbuf);
1856     }
1857     /* These switches take effect only at launch time */
1858     if (FALSE == bAppendFiles)
1859     {
1860         add_to_string(cmd_args_launch, "-noappend ");
1861     }
1862     if (bKeepAndNumCPT)
1863     {
1864         add_to_string(cmd_args_launch, "-cpnum ");
1865     }
1866     if (bResetHWay)
1867     {
1868         add_to_string(cmd_args_launch, "-resethway ");
1869     }
1870
1871     /********************/
1872     /* 2. Process files */
1873     /********************/
1874     for (i = 0; i < nfile; i++)
1875     {
1876         opt  = (char *)fnm[i].opt;
1877         name = opt2fn(opt, nfile, fnm);
1878
1879         /* Strbuf contains the options, now let's sort out where we need that */
1880         sprintf(strbuf, "%s %s ", opt, name);
1881
1882         if (is_bench_file(opt, opt2bSet(opt, nfile, fnm), is_optional(&fnm[i]), is_output(&fnm[i])) )
1883         {
1884             /* All options starting with -b* need the 'b' removed,
1885              * therefore overwrite strbuf */
1886             if (0 == strncmp(opt, "-b", 2))
1887             {
1888                 sprintf(strbuf, "-%s %s ", &opt[2], name);
1889             }
1890
1891             add_to_string(cmd_args_bench, strbuf);
1892         }
1893
1894         if (is_launch_file(opt, opt2bSet(opt, nfile, fnm)) )
1895         {
1896             add_to_string(cmd_args_launch, strbuf);
1897         }
1898     }
1899
1900     add_to_string(cmd_args_bench, extra_args);
1901     add_to_string(cmd_args_launch, extra_args);
1902 }
1903
1904
1905 /* Set option opt */
1906 static void setopt(const char *opt, int nfile, t_filenm fnm[])
1907 {
1908     int i;
1909
1910     for (i = 0; (i < nfile); i++)
1911     {
1912         if (strcmp(opt, fnm[i].opt) == 0)
1913         {
1914             fnm[i].flag |= ffSET;
1915         }
1916     }
1917 }
1918
1919
1920 /* This routine inspects the tpr file and ...
1921  * 1. checks for output files that get triggered by a tpr option. These output
1922  *    files are marked as 'set' to allow for a proper cleanup after each
1923  *    tuning run.
1924  * 2. returns the PME:PP load ratio
1925  * 3. returns rcoulomb from the tpr */
1926 static float inspect_tpr(int nfile, t_filenm fnm[], real *rcoulomb)
1927 {
1928     gmx_bool     bPull;     /* Is pulling requested in .tpr file?             */
1929     gmx_bool     bTpi;      /* Is test particle insertion requested?          */
1930     gmx_bool     bFree;     /* Is a free energy simulation requested?         */
1931     gmx_bool     bNM;       /* Is a normal mode analysis requested?           */
1932     gmx_bool     bSwap;     /* Is water/ion position swapping requested?      */
1933     t_inputrec   ir;
1934     t_state      state;
1935     gmx_mtop_t   mtop;
1936
1937
1938     /* Check tpr file for options that trigger extra output files */
1939     read_tpx_state(opt2fn("-s", nfile, fnm), &ir, &state, NULL, &mtop);
1940     bPull = (epullNO != ir.ePull);
1941     bFree = (efepNO  != ir.efep );
1942     bNM   = (eiNM    == ir.eI   );
1943     bSwap = (eswapNO != ir.eSwapCoords);
1944     bTpi  = EI_TPI(ir.eI);
1945
1946     /* Set these output files on the tuning command-line */
1947     if (bPull)
1948     {
1949         setopt("-pf", nfile, fnm);
1950         setopt("-px", nfile, fnm);
1951     }
1952     if (bFree)
1953     {
1954         setopt("-dhdl", nfile, fnm);
1955     }
1956     if (bTpi)
1957     {
1958         setopt("-tpi", nfile, fnm);
1959         setopt("-tpid", nfile, fnm);
1960     }
1961     if (bNM)
1962     {
1963         setopt("-mtx", nfile, fnm);
1964     }
1965     if (bSwap)
1966     {
1967         setopt("-swap", nfile, fnm);
1968     }
1969
1970     *rcoulomb = ir.rcoulomb;
1971
1972     /* Return the estimate for the number of PME nodes */
1973     return pme_load_estimate(&mtop, &ir, state.box);
1974 }
1975
1976
1977 static void couple_files_options(int nfile, t_filenm fnm[])
1978 {
1979     int      i;
1980     gmx_bool bSet, bBench;
1981     char    *opt;
1982     char     buf[20];
1983
1984
1985     for (i = 0; i < nfile; i++)
1986     {
1987         opt    = (char *)fnm[i].opt;
1988         bSet   = ((fnm[i].flag & ffSET) != 0);
1989         bBench = (0 == strncmp(opt, "-b", 2));
1990
1991         /* Check optional files */
1992         /* If e.g. -eo is set, then -beo also needs to be set */
1993         if (is_optional(&fnm[i]) && bSet && !bBench)
1994         {
1995             sprintf(buf, "-b%s", &opt[1]);
1996             setopt(buf, nfile, fnm);
1997         }
1998         /* If -beo is set, then -eo also needs to be! */
1999         if (is_optional(&fnm[i]) && bSet && bBench)
2000         {
2001             sprintf(buf, "-%s", &opt[2]);
2002             setopt(buf, nfile, fnm);
2003         }
2004     }
2005 }
2006
2007
2008 #define BENCHSTEPS (1000)
2009
2010 int gmx_tune_pme(int argc, char *argv[])
2011 {
2012     const char     *desc[] = {
2013         "For a given number [TT]-np[tt] or [TT]-ntmpi[tt] of ranks, [THISMODULE] systematically",
2014         "times [gmx-mdrun] with various numbers of PME-only ranks and determines",
2015         "which setting is fastest. It will also test whether performance can",
2016         "be enhanced by shifting load from the reciprocal to the real space",
2017         "part of the Ewald sum. ",
2018         "Simply pass your [TT].tpr[tt] file to [THISMODULE] together with other options",
2019         "for [gmx-mdrun] as needed.[PAR]",
2020         "Which executables are used can be set in the environment variables",
2021         "MPIRUN and MDRUN. If these are not present, 'mpirun' and 'mdrun'",
2022         "will be used as defaults. Note that for certain MPI frameworks you",
2023         "need to provide a machine- or hostfile. This can also be passed",
2024         "via the MPIRUN variable, e.g.[PAR]",
2025         "[TT]export MPIRUN=\"/usr/local/mpirun -machinefile hosts\"[tt][PAR]",
2026         "Before doing the actual benchmark runs, [THISMODULE] will do a quick",
2027         "check whether mdrun works as expected with the provided parallel settings",
2028         "if the [TT]-check[tt] option is activated (the default).",
2029         "Please call [THISMODULE] with the normal options you would pass to",
2030         "[gmx-mdrun] and add [TT]-np[tt] for the number of ranks to perform the",
2031         "tests on, or [TT]-ntmpi[tt] for the number of threads. You can also add [TT]-r[tt]",
2032         "to repeat each test several times to get better statistics. [PAR]",
2033         "[THISMODULE] can test various real space / reciprocal space workloads",
2034         "for you. With [TT]-ntpr[tt] you control how many extra [TT].tpr[tt] files will be",
2035         "written with enlarged cutoffs and smaller Fourier grids respectively.",
2036         "Typically, the first test (number 0) will be with the settings from the input",
2037         "[TT].tpr[tt] file; the last test (number [TT]ntpr[tt]) will have the Coulomb cutoff",
2038         "specified by [TT]-rmax[tt] with a somewhat smaller PME grid at the same time. ",
2039         "In this last test, the Fourier spacing is multiplied with [TT]rmax[tt]/rcoulomb. ",
2040         "The remaining [TT].tpr[tt] files will have equally-spaced Coulomb radii (and Fourier "
2041         "spacings) between these extremes. [BB]Note[bb] that you can set [TT]-ntpr[tt] to 1",
2042         "if you just seek the optimal number of PME-only ranks; in that case",
2043         "your input [TT].tpr[tt] file will remain unchanged.[PAR]",
2044         "For the benchmark runs, the default of 1000 time steps should suffice for most",
2045         "MD systems. The dynamic load balancing needs about 100 time steps",
2046         "to adapt to local load imbalances, therefore the time step counters",
2047         "are by default reset after 100 steps. For large systems (>1M atoms), as well as ",
2048         "for a higher accuracy of the measurements, you should set [TT]-resetstep[tt] to a higher value.",
2049         "From the 'DD' load imbalance entries in the md.log output file you",
2050         "can tell after how many steps the load is sufficiently balanced. Example call:[PAR]"
2051         "[TT]gmx tune_pme -np 64 -s protein.tpr -launch[tt][PAR]",
2052         "After calling [gmx-mdrun] several times, detailed performance information",
2053         "is available in the output file [TT]perf.out[tt].",
2054         "[BB]Note[bb] that during the benchmarks, a couple of temporary files are written",
2055         "(options [TT]-b[tt]*), these will be automatically deleted after each test.[PAR]",
2056         "If you want the simulation to be started automatically with the",
2057         "optimized parameters, use the command line option [TT]-launch[tt].[PAR]",
2058     };
2059
2060     int             nnodes         = 1;
2061     int             repeats        = 2;
2062     int             pmeentries     = 0; /* How many values for -npme do we actually test for each tpr file */
2063     real            maxPMEfraction = 0.50;
2064     real            minPMEfraction = 0.25;
2065     int             maxPMEnodes, minPMEnodes;
2066     float           guessPMEratio;                    /* guessed PME:PP ratio based on the tpr file */
2067     float           guessPMEnodes;
2068     int             npme_fixed     = -2;              /* If >= -1, use only this number
2069                                                        * of PME-only nodes                */
2070     int             ntprs          = 0;
2071     real            rmin           = 0.0, rmax = 0.0; /* min and max value for rcoulomb if scaling is requested */
2072     real            rcoulomb       = -1.0;            /* Coulomb radius as set in .tpr file */
2073     gmx_bool        bScaleRvdw     = TRUE;
2074     gmx_int64_t     bench_nsteps   = BENCHSTEPS;
2075     gmx_int64_t     new_sim_nsteps = -1;  /* -1 indicates: not set by the user */
2076     gmx_int64_t     cpt_steps      = 0;   /* Step counter in .cpt input file   */
2077     int             presteps       = 100; /* Do a full cycle reset after presteps steps */
2078     gmx_bool        bOverwrite     = FALSE, bKeepTPR;
2079     gmx_bool        bLaunch        = FALSE;
2080     char           *ExtraArgs      = NULL;
2081     char          **tpr_names      = NULL;
2082     const char     *simulation_tpr = NULL;
2083     int             best_npme, best_tpr;
2084     int             sim_part = 1; /* For benchmarks with checkpoint files */
2085     char            bbuf[STRLEN];
2086
2087     /* Default program names if nothing else is found */
2088     char         *cmd_mpirun = NULL, *cmd_mdrun = NULL;
2089     char         *cmd_args_bench, *cmd_args_launch;
2090     char         *cmd_np = NULL;
2091
2092     t_perf      **perfdata = NULL;
2093     t_inputinfo  *info;
2094     int           i;
2095     FILE         *fp;
2096     t_commrec    *cr;
2097
2098     /* Print out how long the tuning took */
2099     double          seconds;
2100
2101     static t_filenm fnm[] = {
2102         /* g_tune_pme */
2103         { efOUT, "-p",      "perf",     ffWRITE },
2104         { efLOG, "-err",    "bencherr", ffWRITE },
2105         { efTPR, "-so",     "tuned",    ffWRITE },
2106         /* mdrun: */
2107         { efTPR, NULL,      NULL,       ffREAD },
2108         { efTRN, "-o",      NULL,       ffWRITE },
2109         { efCOMPRESSED, "-x", NULL,     ffOPTWR },
2110         { efCPT, "-cpi",    NULL,       ffOPTRD },
2111         { efCPT, "-cpo",    NULL,       ffOPTWR },
2112         { efSTO, "-c",      "confout",  ffWRITE },
2113         { efEDR, "-e",      "ener",     ffWRITE },
2114         { efLOG, "-g",      "md",       ffWRITE },
2115         { efXVG, "-dhdl",   "dhdl",     ffOPTWR },
2116         { efXVG, "-field",  "field",    ffOPTWR },
2117         { efXVG, "-table",  "table",    ffOPTRD },
2118         { efXVG, "-tabletf", "tabletf",   ffOPTRD },
2119         { efXVG, "-tablep", "tablep",   ffOPTRD },
2120         { efXVG, "-tableb", "table",    ffOPTRD },
2121         { efTRX, "-rerun",  "rerun",    ffOPTRD },
2122         { efXVG, "-tpi",    "tpi",      ffOPTWR },
2123         { efXVG, "-tpid",   "tpidist",  ffOPTWR },
2124         { efEDI, "-ei",     "sam",      ffOPTRD },
2125         { efXVG, "-eo",     "edsam",    ffOPTWR },
2126         { efXVG, "-devout", "deviatie", ffOPTWR },
2127         { efXVG, "-runav",  "runaver",  ffOPTWR },
2128         { efXVG, "-px",     "pullx",    ffOPTWR },
2129         { efXVG, "-pf",     "pullf",    ffOPTWR },
2130         { efXVG, "-ro",     "rotation", ffOPTWR },
2131         { efLOG, "-ra",     "rotangles", ffOPTWR },
2132         { efLOG, "-rs",     "rotslabs", ffOPTWR },
2133         { efLOG, "-rt",     "rottorque", ffOPTWR },
2134         { efMTX, "-mtx",    "nm",       ffOPTWR },
2135         { efNDX, "-dn",     "dipole",   ffOPTWR },
2136         { efXVG, "-swap",   "swapions", ffOPTWR },
2137         /* Output files that are deleted after each benchmark run */
2138         { efTRN, "-bo",     "bench",    ffWRITE },
2139         { efXTC, "-bx",     "bench",    ffWRITE },
2140         { efCPT, "-bcpo",   "bench",    ffWRITE },
2141         { efSTO, "-bc",     "bench",    ffWRITE },
2142         { efEDR, "-be",     "bench",    ffWRITE },
2143         { efLOG, "-bg",     "bench",    ffWRITE },
2144         { efXVG, "-beo",    "benchedo", ffOPTWR },
2145         { efXVG, "-bdhdl",  "benchdhdl", ffOPTWR },
2146         { efXVG, "-bfield", "benchfld", ffOPTWR },
2147         { efXVG, "-btpi",   "benchtpi", ffOPTWR },
2148         { efXVG, "-btpid",  "benchtpid", ffOPTWR },
2149         { efXVG, "-bdevout", "benchdev", ffOPTWR },
2150         { efXVG, "-brunav", "benchrnav", ffOPTWR },
2151         { efXVG, "-bpx",    "benchpx",  ffOPTWR },
2152         { efXVG, "-bpf",    "benchpf",  ffOPTWR },
2153         { efXVG, "-bro",    "benchrot", ffOPTWR },
2154         { efLOG, "-bra",    "benchrota", ffOPTWR },
2155         { efLOG, "-brs",    "benchrots", ffOPTWR },
2156         { efLOG, "-brt",    "benchrott", ffOPTWR },
2157         { efMTX, "-bmtx",   "benchn",   ffOPTWR },
2158         { efNDX, "-bdn",    "bench",    ffOPTWR },
2159         { efXVG, "-bswap",  "benchswp", ffOPTWR }
2160     };
2161
2162     gmx_bool        bThreads     = FALSE;
2163
2164     int             nthreads = 1;
2165
2166     const char     *procstring[] =
2167     { NULL, "-np", "-n", "none", NULL };
2168     const char     *npmevalues_opt[] =
2169     { NULL, "auto", "all", "subset", NULL };
2170
2171     gmx_bool     bAppendFiles          = TRUE;
2172     gmx_bool     bKeepAndNumCPT        = FALSE;
2173     gmx_bool     bResetCountersHalfWay = FALSE;
2174     gmx_bool     bBenchmark            = TRUE;
2175     gmx_bool     bCheck                = TRUE;
2176
2177     output_env_t oenv = NULL;
2178
2179     t_pargs      pa[] = {
2180         /***********************/
2181         /* g_tune_pme options: */
2182         /***********************/
2183         { "-np",       FALSE, etINT,  {&nnodes},
2184           "Number of ranks to run the tests on (must be > 2 for separate PME ranks)" },
2185         { "-npstring", FALSE, etENUM, {procstring},
2186           "Specify the number of ranks to [TT]$MPIRUN[tt] using this string"},
2187         { "-ntmpi",    FALSE, etINT,  {&nthreads},
2188           "Number of MPI-threads to run the tests on (turns MPI & mpirun off)"},
2189         { "-r",        FALSE, etINT,  {&repeats},
2190           "Repeat each test this often" },
2191         { "-max",      FALSE, etREAL, {&maxPMEfraction},
2192           "Max fraction of PME ranks to test with" },
2193         { "-min",      FALSE, etREAL, {&minPMEfraction},
2194           "Min fraction of PME ranks to test with" },
2195         { "-npme",     FALSE, etENUM, {npmevalues_opt},
2196           "Within -min and -max, benchmark all possible values for [TT]-npme[tt], or just a reasonable subset. "
2197           "Auto neglects -min and -max and chooses reasonable values around a guess for npme derived from the .tpr"},
2198         { "-fix",      FALSE, etINT,  {&npme_fixed},
2199           "If >= -1, do not vary the number of PME-only ranks, instead use this fixed value and only vary rcoulomb and the PME grid spacing."},
2200         { "-rmax",     FALSE, etREAL, {&rmax},
2201           "If >0, maximal rcoulomb for -ntpr>1 (rcoulomb upscaling results in fourier grid downscaling)" },
2202         { "-rmin",     FALSE, etREAL, {&rmin},
2203           "If >0, minimal rcoulomb for -ntpr>1" },
2204         { "-scalevdw",  FALSE, etBOOL, {&bScaleRvdw},
2205           "Scale rvdw along with rcoulomb"},
2206         { "-ntpr",     FALSE, etINT,  {&ntprs},
2207           "Number of [TT].tpr[tt] files to benchmark. Create this many files with different rcoulomb scaling factors depending on -rmin and -rmax. "
2208           "If < 1, automatically choose the number of [TT].tpr[tt] files to test" },
2209         { "-steps",    FALSE, etINT64, {&bench_nsteps},
2210           "Take timings for this many steps in the benchmark runs" },
2211         { "-resetstep", FALSE, etINT,  {&presteps},
2212           "Let dlb equilibrate this many steps before timings are taken (reset cycle counters after this many steps)" },
2213         { "-simsteps", FALSE, etINT64, {&new_sim_nsteps},
2214           "If non-negative, perform this many steps in the real run (overwrites nsteps from [TT].tpr[tt], add [TT].cpt[tt] steps)" },
2215         { "-launch",   FALSE, etBOOL, {&bLaunch},
2216           "Launch the real simulation after optimization" },
2217         { "-bench",    FALSE, etBOOL, {&bBenchmark},
2218           "Run the benchmarks or just create the input [TT].tpr[tt] files?" },
2219         { "-check",    FALSE, etBOOL, {&bCheck},
2220           "Before the benchmark runs, check whether mdrun works in parallel" },
2221         /******************/
2222         /* mdrun options: */
2223         /******************/
2224         /* We let g_tune_pme parse and understand these options, because we need to
2225          * prevent that they appear on the mdrun command line for the benchmarks */
2226         { "-append",   FALSE, etBOOL, {&bAppendFiles},
2227           "Append to previous output files when continuing from checkpoint instead of adding the simulation part number to all file names (for launch only)" },
2228         { "-cpnum",    FALSE, etBOOL, {&bKeepAndNumCPT},
2229           "Keep and number checkpoint files (launch only)" },
2230         { "-resethway", FALSE, etBOOL, {&bResetCountersHalfWay},
2231           "HIDDENReset the cycle counters after half the number of steps or halfway [TT]-maxh[tt] (launch only)" }
2232     };
2233
2234 #define NFILE asize(fnm)
2235
2236     seconds = gmx_gettime();
2237
2238     if (!parse_common_args(&argc, argv, PCA_NOEXIT_ON_ARGS,
2239                            NFILE, fnm, asize(pa), pa, asize(desc), desc,
2240                            0, NULL, &oenv))
2241     {
2242         return 0;
2243     }
2244
2245     /* Store the remaining unparsed command line entries in a string which
2246      * is then attached to the mdrun command line */
2247     snew(ExtraArgs, 1);
2248     ExtraArgs[0] = '\0';
2249     for (i = 1; i < argc; i++) /* argc will now be 1 if everything was understood */
2250     {
2251         add_to_string(&ExtraArgs, argv[i]);
2252         add_to_string(&ExtraArgs, " ");
2253     }
2254
2255     if (opt2parg_bSet("-ntmpi", asize(pa), pa))
2256     {
2257         bThreads = TRUE;
2258         if (opt2parg_bSet("-npstring", asize(pa), pa))
2259         {
2260             fprintf(stderr, "WARNING: -npstring has no effect when using threads.\n");
2261         }
2262
2263         if (nnodes > 1)
2264         {
2265             gmx_fatal(FARGS, "Can't run multi-threaded MPI simulation yet!");
2266         }
2267         /* and now we just set this; a bit of an ugly hack*/
2268         nnodes = nthreads;
2269     }
2270     /* Check for PME:PP ratio and whether tpr triggers additional output files */
2271     guessPMEratio = inspect_tpr(NFILE, fnm, &rcoulomb);
2272
2273     /* Automatically set -beo options if -eo is set etc. */
2274     couple_files_options(NFILE, fnm);
2275
2276     /* Construct the command line arguments for benchmark runs
2277      * as well as for the simulation run */
2278     if (bThreads)
2279     {
2280         sprintf(bbuf, " -ntmpi %d ", nthreads);
2281     }
2282     else
2283     {
2284         /* This string will be used for MPI runs and will appear after the
2285          * mpirun command. */
2286         if (strcmp(procstring[0], "none") != 0)
2287         {
2288             sprintf(bbuf, " %s %d ", procstring[0], nnodes);
2289         }
2290         else
2291         {
2292             sprintf(bbuf, " ");
2293         }
2294     }
2295
2296     cmd_np = bbuf;
2297
2298     create_command_line_snippets(bAppendFiles, bKeepAndNumCPT, bResetCountersHalfWay, presteps,
2299                                  NFILE, fnm, &cmd_args_bench, &cmd_args_launch, ExtraArgs);
2300
2301     /* Read in checkpoint file if requested */
2302     sim_part = 1;
2303     if (opt2bSet("-cpi", NFILE, fnm))
2304     {
2305         snew(cr, 1);
2306         cr->duty = DUTY_PP; /* makes the following routine happy */
2307         read_checkpoint_simulation_part(opt2fn("-cpi", NFILE, fnm),
2308                                         &sim_part, &cpt_steps, cr,
2309                                         FALSE, NFILE, fnm, NULL, NULL);
2310         sfree(cr);
2311         sim_part++;
2312         /* sim_part will now be 1 if no checkpoint file was found */
2313         if (sim_part <= 1)
2314         {
2315             gmx_fatal(FARGS, "Checkpoint file %s not found!", opt2fn("-cpi", NFILE, fnm));
2316         }
2317     }
2318
2319     /* Open performance output file and write header info */
2320     fp = gmx_ffopen(opt2fn("-p", NFILE, fnm), "w");
2321
2322     /* Make a quick consistency check of command line parameters */
2323     check_input(nnodes, repeats, &ntprs, &rmin, rcoulomb, &rmax,
2324                 maxPMEfraction, minPMEfraction, npme_fixed,
2325                 bench_nsteps, fnm, NFILE, sim_part, presteps,
2326                 asize(pa), pa);
2327
2328     /* Determine the maximum and minimum number of PME nodes to test,
2329      * the actual list of settings is build in do_the_tests(). */
2330     if ((nnodes > 2) && (npme_fixed < -1))
2331     {
2332         if (0 == strcmp(npmevalues_opt[0], "auto"))
2333         {
2334             /* Determine the npme range automatically based on the PME:PP load guess */
2335             if (guessPMEratio > 1.0)
2336             {
2337                 /* More PME than PP work, probably we do not need separate PME nodes at all! */
2338                 maxPMEnodes = nnodes/2;
2339                 minPMEnodes = nnodes/2;
2340             }
2341             else
2342             {
2343                 /* PME : PP load is in the range 0..1, let's test around the guess */
2344                 guessPMEnodes = nnodes/(1.0 + 1.0/guessPMEratio);
2345                 minPMEnodes   = floor(0.7*guessPMEnodes);
2346                 maxPMEnodes   =  ceil(1.6*guessPMEnodes);
2347                 maxPMEnodes   = min(maxPMEnodes, nnodes/2);
2348             }
2349         }
2350         else
2351         {
2352             /* Determine the npme range based on user input */
2353             maxPMEnodes = floor(maxPMEfraction*nnodes);
2354             minPMEnodes = max(floor(minPMEfraction*nnodes), 0);
2355             fprintf(stdout, "Will try runs with %d ", minPMEnodes);
2356             if (maxPMEnodes != minPMEnodes)
2357             {
2358                 fprintf(stdout, "- %d ", maxPMEnodes);
2359             }
2360             fprintf(stdout, "PME-only ranks.\n  Note that the automatic number of PME-only ranks and no separate PME ranks are always tested.\n");
2361         }
2362     }
2363     else
2364     {
2365         maxPMEnodes = 0;
2366         minPMEnodes = 0;
2367     }
2368
2369     /* Get the commands we need to set up the runs from environment variables */
2370     get_program_paths(bThreads, &cmd_mpirun, &cmd_mdrun);
2371     if (bBenchmark && repeats > 0)
2372     {
2373         check_mdrun_works(bThreads, cmd_mpirun, cmd_np, cmd_mdrun);
2374     }
2375
2376     /* Print some header info to file */
2377     sep_line(fp);
2378     fprintf(fp, "\n      P E R F O R M A N C E   R E S U L T S\n");
2379     sep_line(fp);
2380     fprintf(fp, "%s for Gromacs %s\n", output_env_get_program_display_name(oenv),
2381             gmx_version());
2382     if (!bThreads)
2383     {
2384         fprintf(fp, "Number of ranks         : %d\n", nnodes);
2385         fprintf(fp, "The mpirun command is   : %s\n", cmd_mpirun);
2386         if (strcmp(procstring[0], "none") != 0)
2387         {
2388             fprintf(fp, "Passing # of ranks via  : %s\n", procstring[0]);
2389         }
2390         else
2391         {
2392             fprintf(fp, "Not setting number of ranks in system call\n");
2393         }
2394     }
2395     else
2396     {
2397         fprintf(fp, "Number of threads       : %d\n", nnodes);
2398     }
2399
2400     fprintf(fp, "The mdrun  command is   : %s\n", cmd_mdrun);
2401     fprintf(fp, "mdrun args benchmarks   : %s\n", cmd_args_bench);
2402     fprintf(fp, "Benchmark steps         : ");
2403     fprintf(fp, "%"GMX_PRId64, bench_nsteps);
2404     fprintf(fp, "\n");
2405     fprintf(fp, "dlb equilibration steps : %d\n", presteps);
2406     if (sim_part > 1)
2407     {
2408         fprintf(fp, "Checkpoint time step    : ");
2409         fprintf(fp, "%"GMX_PRId64, cpt_steps);
2410         fprintf(fp, "\n");
2411     }
2412     fprintf(fp, "mdrun args at launchtime: %s\n", cmd_args_launch);
2413
2414     if (new_sim_nsteps >= 0)
2415     {
2416         bOverwrite = TRUE;
2417         fprintf(stderr, "Note: Simulation input file %s will have ", opt2fn("-so", NFILE, fnm));
2418         fprintf(stderr, "%"GMX_PRId64, new_sim_nsteps+cpt_steps);
2419         fprintf(stderr, " steps.\n");
2420         fprintf(fp, "Simulation steps        : ");
2421         fprintf(fp, "%"GMX_PRId64, new_sim_nsteps);
2422         fprintf(fp, "\n");
2423     }
2424     if (repeats > 1)
2425     {
2426         fprintf(fp, "Repeats for each test   : %d\n", repeats);
2427     }
2428
2429     if (npme_fixed >= -1)
2430     {
2431         fprintf(fp, "Fixing -npme at         : %d\n", npme_fixed);
2432     }
2433
2434     fprintf(fp, "Input file              : %s\n", opt2fn("-s", NFILE, fnm));
2435     fprintf(fp, "   PME/PP load estimate : %g\n", guessPMEratio);
2436
2437     /* Allocate memory for the inputinfo struct: */
2438     snew(info, 1);
2439     info->nr_inputfiles = ntprs;
2440     for (i = 0; i < ntprs; i++)
2441     {
2442         snew(info->rcoulomb, ntprs);
2443         snew(info->rvdw, ntprs);
2444         snew(info->rlist, ntprs);
2445         snew(info->rlistlong, ntprs);
2446         snew(info->nkx, ntprs);
2447         snew(info->nky, ntprs);
2448         snew(info->nkz, ntprs);
2449         snew(info->fsx, ntprs);
2450         snew(info->fsy, ntprs);
2451         snew(info->fsz, ntprs);
2452     }
2453     /* Make alternative tpr files to test: */
2454     snew(tpr_names, ntprs);
2455     for (i = 0; i < ntprs; i++)
2456     {
2457         snew(tpr_names[i], STRLEN);
2458     }
2459
2460     /* It can be that ntprs is reduced by make_benchmark_tprs if not enough
2461      * different grids could be found. */
2462     make_benchmark_tprs(opt2fn("-s", NFILE, fnm), tpr_names, bench_nsteps+presteps,
2463                         cpt_steps, rmin, rmax, bScaleRvdw, &ntprs, info, fp);
2464
2465     /********************************************************************************/
2466     /* Main loop over all scenarios we need to test: tpr files, PME nodes, repeats  */
2467     /********************************************************************************/
2468     snew(perfdata, ntprs);
2469     if (bBenchmark)
2470     {
2471         do_the_tests(fp, tpr_names, maxPMEnodes, minPMEnodes, npme_fixed, npmevalues_opt[0], perfdata, &pmeentries,
2472                      repeats, nnodes, ntprs, bThreads, cmd_mpirun, cmd_np, cmd_mdrun,
2473                      cmd_args_bench, fnm, NFILE, presteps, cpt_steps, bCheck);
2474
2475         fprintf(fp, "\nTuning took%8.1f minutes.\n", (gmx_gettime()-seconds)/60.0);
2476
2477         /* Analyse the results and give a suggestion for optimal settings: */
2478         bKeepTPR = analyze_data(fp, opt2fn("-p", NFILE, fnm), perfdata, nnodes, ntprs, pmeentries,
2479                                 repeats, info, &best_tpr, &best_npme);
2480
2481         /* Take the best-performing tpr file and enlarge nsteps to original value */
2482         if (bKeepTPR && !bOverwrite)
2483         {
2484             simulation_tpr = opt2fn("-s", NFILE, fnm);
2485         }
2486         else
2487         {
2488             simulation_tpr = opt2fn("-so", NFILE, fnm);
2489             modify_PMEsettings(bOverwrite ? (new_sim_nsteps+cpt_steps) : info->orig_sim_steps,
2490                                info->orig_init_step, tpr_names[best_tpr], simulation_tpr);
2491         }
2492
2493         /* Let's get rid of the temporary benchmark input files */
2494         for (i = 0; i < ntprs; i++)
2495         {
2496             fprintf(stdout, "Deleting temporary benchmark input file %s\n", tpr_names[i]);
2497             remove(tpr_names[i]);
2498         }
2499
2500         /* Now start the real simulation if the user requested it ... */
2501         launch_simulation(bLaunch, fp, bThreads, cmd_mpirun, cmd_np, cmd_mdrun,
2502                           cmd_args_launch, simulation_tpr, best_npme);
2503     }
2504     gmx_ffclose(fp);
2505
2506     /* ... or simply print the performance results to screen: */
2507     if (!bLaunch)
2508     {
2509         finalize(opt2fn("-p", NFILE, fnm));
2510     }
2511
2512     return 0;
2513 }