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