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