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