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