src/tools/gmx_tune_pme.c

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