src/tools/gmx_tune_pme.c

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