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