2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016,2017,2018,2019, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
47 #include "gromacs/commandline/pargs.h"
48 #include "gromacs/fileio/enxio.h"
49 #include "gromacs/fileio/trxio.h"
50 #include "gromacs/listed_forces/disre.h"
51 #include "gromacs/math/functions.h"
52 #include "gromacs/math/vec.h"
53 #include "gromacs/mdtypes/md_enums.h"
54 #include "gromacs/trajectory/energyframe.h"
55 #include "gromacs/utility/arrayref.h"
56 #include "gromacs/utility/arraysize.h"
57 #include "gromacs/utility/cstringutil.h"
58 #include "gromacs/utility/fatalerror.h"
59 #include "gromacs/utility/int64_to_int.h"
60 #include "gromacs/utility/smalloc.h"
61 #include "gromacs/utility/strconvert.h"
63 #define TIME_EXPLICIT 0
64 #define TIME_CONTINUE 1
70 static int* select_it(int nre, gmx_enxnm_t* nm, int* nset)
75 gmx_bool bVerbose = TRUE;
77 if ((getenv("GMX_ENER_VERBOSE")) != nullptr)
82 fprintf(stderr, "Select the terms you want to scale from the following list\n");
83 fprintf(stderr, "End your selection with 0\n");
87 for (k = 0; (k < nre);)
89 for (j = 0; (j < 4) && (k < nre); j++, k++)
91 fprintf(stderr, " %3d=%14s", k + 1, nm[k].name);
93 fprintf(stderr, "\n");
100 if (1 != scanf("%d", &n))
102 gmx_fatal(FARGS, "Cannot read energy term");
104 if ((n > 0) && (n <= nre))
111 for (i = (*nset) = 0; (i < nre); i++)
124 static void sort_files(gmx::ArrayRef<std::string> files, real* settime)
126 for (gmx::index i = 0; i < files.ssize(); i++)
128 gmx::index minidx = i;
129 for (gmx::index j = i + 1; j < files.ssize(); j++)
131 if (settime[j] < settime[minidx])
138 real timeswap = settime[i];
139 settime[i] = settime[minidx];
140 settime[minidx] = timeswap;
141 std::string tmp = files[i];
142 files[i] = files[minidx];
149 static int scan_ene_files(const std::vector<std::string>& files, real* readtime, real* timestep, int* nremax)
151 /* Check number of energy terms and start time of all files */
152 int nre, nremin = 0, nresav = 0;
155 char inputstring[STRLEN];
161 for (size_t f = 0; f < files.size(); f++)
163 in = open_enx(files[f].c_str(), "r");
165 do_enxnms(in, &nre, &enm);
182 nremin = std::min(nremin, fr->nre);
183 *nremax = std::max(*nremax, fr->nre);
187 "Energy files don't match, different number of energies:\n"
188 " %s: %d\n %s: %d\n",
189 files[f - 1].c_str(), nresav, files[f].c_str(), fr->nre);
191 "\nContinue conversion using only the first %d terms (n/y)?\n"
192 "(you should be sure that the energy terms match)\n",
194 if (nullptr == fgets(inputstring, STRLEN - 1, stdin))
196 gmx_fatal(FARGS, "Error reading user input");
198 if (inputstring[0] != 'y' && inputstring[0] != 'Y')
200 fprintf(stderr, "Will not convert\n");
209 fprintf(stderr, "\n");
210 free_enxnms(nre, enm);
220 static void edit_files(gmx::ArrayRef<std::string> files,
228 char inputstring[STRLEN], *chptr;
232 if (files.size() == 1)
234 fprintf(stderr, "\n\nEnter the new start time:\n\n");
239 "\n\nEnter the new start time for each file.\n"
240 "There are two special options, both disables sorting:\n\n"
241 "c (continue) - The start time is taken from the end\n"
242 "of the previous file. Use it when your continuation run\n"
243 "restarts with t=0 and there is no overlap.\n\n"
244 "l (last) - The time in this file will be changed the\n"
245 "same amount as in the previous. Use it when the time in the\n"
246 "new run continues from the end of the previous one,\n"
247 "since this takes possible overlap into account.\n\n");
251 " File Current start New start\n"
252 "---------------------------------------------------------\n");
254 for (gmx::index i = 0; i < files.ssize(); i++)
256 fprintf(stderr, "%25s %10.3f ", files[i].c_str(), readtime[i]);
260 if (nullptr == fgets(inputstring, STRLEN - 1, stdin))
262 gmx_fatal(FARGS, "Error reading user input");
264 inputstring[std::strlen(inputstring) - 1] = 0;
266 if (inputstring[0] == 'c' || inputstring[0] == 'C')
268 cont_type[i] = TIME_CONTINUE;
271 settime[i] = FLT_MAX;
273 else if (inputstring[0] == 'l' || inputstring[0] == 'L')
275 cont_type[i] = TIME_LAST;
278 settime[i] = FLT_MAX;
282 settime[i] = strtod(inputstring, &chptr);
283 if (chptr == inputstring)
285 fprintf(stderr, "Try that again: ");
289 cont_type[i] = TIME_EXPLICIT;
295 if (cont_type[0] != TIME_EXPLICIT)
297 cont_type[0] = TIME_EXPLICIT;
303 for (gmx::index i = 0; i < files.ssize(); i++)
305 settime[i] = readtime[i];
309 if (bSort && files.size() > 1)
311 sort_files(files, settime);
315 fprintf(stderr, "Sorting disabled.\n");
319 /* Write out the new order and start times */
321 "\nSummary of files and start times used:\n\n"
323 "-----------------------------------------\n");
324 for (gmx::index i = 0; i < files.ssize(); i++)
326 switch (cont_type[i])
329 fprintf(stderr, "%25s %10.3f\n", files[i].c_str(), settime[i]);
332 fprintf(stderr, "%25s Continue from end of last file\n", files[i].c_str());
335 fprintf(stderr, "%25s Change by same amount as last file\n", files[i].c_str());
339 fprintf(stderr, "\n");
341 settime[files.size()] = FLT_MAX;
342 cont_type[files.size()] = TIME_EXPLICIT;
343 readtime[files.size()] = FLT_MAX;
347 static void update_ee_sum(int nre,
348 int64_t* ee_sum_step,
349 int64_t* ee_sum_nsteps,
350 int64_t* ee_sum_nsum,
355 int64_t nsteps, nsum, fr_nsum;
358 nsteps = *ee_sum_nsteps;
371 for (i = 0; i < nre; i++)
373 ee_sum[i].esum = fr->ener[i].e;
379 for (i = 0; i < nre; i++)
381 ee_sum[i].esum = fr->ener[i].esum;
382 ee_sum[i].eav = fr->ener[i].eav;
388 else if (out_step + *ee_sum_nsum - *ee_sum_step == nsteps + fr->nsteps)
392 for (i = 0; i < nre; i++)
394 ee_sum[i].eav += gmx::square(ee_sum[i].esum / nsum
395 - (ee_sum[i].esum + fr->ener[i].e) / (nsum + 1))
397 ee_sum[i].esum += fr->ener[i].e;
402 for (i = 0; i < fr->nre; i++)
404 ee_sum[i].eav += fr->ener[i].eav
405 + gmx::square(ee_sum[i].esum / nsum
406 - (ee_sum[i].esum + fr->ener[i].esum) / (nsum + fr->nsum))
407 * nsum * (nsum + fr->nsum) / static_cast<double>(fr->nsum);
408 ee_sum[i].esum += fr->ener[i].esum;
411 nsteps += fr->nsteps;
418 fprintf(stderr, "\nWARNING: missing energy sums at time %f\n", fr->t);
424 *ee_sum_step = out_step;
425 *ee_sum_nsteps = nsteps;
429 int gmx_eneconv(int argc, char* argv[])
431 const char* desc[] = {
432 "With [IT]multiple files[it] specified for the [TT]-f[tt] option:[PAR]",
433 "Concatenates several energy files in sorted order.",
434 "In the case of double time frames, the one",
435 "in the later file is used. By specifying [TT]-settime[tt] you will be",
436 "asked for the start time of each file. The input files are taken",
437 "from the command line,",
438 "such that the command [TT]gmx eneconv -f *.edr -o fixed.edr[tt] should do",
440 "With [IT]one file[it] specified for [TT]-f[tt]:[PAR]",
441 "Reads one energy file and writes another, applying the [TT]-dt[tt],",
442 "[TT]-offset[tt], [TT]-t0[tt] and [TT]-settime[tt] options and",
443 "converting to a different format if necessary (indicated by file",
445 "[TT]-settime[tt] is applied first, then [TT]-dt[tt]/[TT]-offset[tt]",
446 "followed by [TT]-b[tt] and [TT]-e[tt] to select which frames to write."
448 const char* bugs[] = {
449 "When combining trajectories the sigma and E^2 (necessary for statistics) are not "
450 "updated correctly. Only the actual energy is correct. One thus has to compute "
451 "statistics in another way."
453 ener_file_t in = nullptr, out = nullptr;
454 gmx_enxnm_t* enm = nullptr;
456 ener_file_t in, out = NULL;
457 gmx_enxnm_t *enm = NULL;
459 t_enxframe * fr, *fro;
460 int64_t ee_sum_step = 0, ee_sum_nsteps, ee_sum_nsum;
462 int64_t lastfilestep, laststep, startstep_file = 0;
464 int nre, nremax, this_nre, i, kkk, nset, *set = nullptr;
466 real * readtime, *settime, timestep, tadjust;
467 char buf[22], buf2[22];
469 gmx_bool bNewFile, bFirst, bNewOutput;
470 gmx_output_env_t* oenv;
471 gmx_bool warned_about_dh = FALSE;
472 t_enxblock* blocks = nullptr;
474 int nblocks_alloc = 0;
477 { efEDR, "-f", nullptr, ffRDMULT },
478 { efEDR, "-o", "fixed", ffWRITE },
481 #define NFILE asize(fnm)
483 static real delta_t = 0.0, toffset = 0, scalefac = 1;
484 static gmx_bool bSetTime = FALSE;
485 static gmx_bool bSort = TRUE, bError = TRUE;
486 static real begin = -1;
487 static real end = -1;
488 gmx_bool remove_dh = FALSE;
491 { "-b", FALSE, etREAL, { &begin }, "First time to use" },
492 { "-e", FALSE, etREAL, { &end }, "Last time to use" },
493 { "-dt", FALSE, etREAL, { &delta_t }, "Only write out frame when t MOD dt = offset" },
494 { "-offset", FALSE, etREAL, { &toffset }, "Time offset for [TT]-dt[tt] option" },
495 { "-settime", FALSE, etBOOL, { &bSetTime }, "Change starting time interactively" },
496 { "-sort", FALSE, etBOOL, { &bSort }, "Sort energy files (not frames)" },
497 { "-rmdh", FALSE, etBOOL, { &remove_dh }, "Remove free energy block data" },
498 { "-scalefac", FALSE, etREAL, { &scalefac }, "Multiply energy component by this factor" },
499 { "-error", FALSE, etBOOL, { &bError }, "Stop on errors in the file" }
502 if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa, asize(desc), desc,
503 asize(bugs), bugs, &oenv))
508 "Note that major changes are planned in future for "
509 "eneconv, to improve usability and utility.");
518 auto files = gmx::copyOf(opt2fns("-f", NFILE, fnm));
522 gmx_fatal(FARGS, "No input files!");
525 snew(settime, files.size() + 1);
526 snew(readtime, files.size() + 1);
527 snew(cont_type, files.size() + 1);
529 nre = scan_ene_files(files, readtime, ×tep, &nremax);
530 edit_files(files, readtime, settime, cont_type, bSetTime, bSort);
534 snew(ee_sum, nremax);
540 snew(fro->ener, nremax);
546 for (size_t f = 0; f < files.size(); f++)
550 in = open_enx(files[f].c_str(), "r");
552 do_enxnms(in, &this_nre, &enm);
557 set = select_it(nre, enm, &nset);
560 /* write names to the output file */
561 out = open_enx(opt2fn("-o", NFILE, fnm), "w");
562 do_enxnms(out, &nre, &enm);
565 /* start reading from the next file */
566 while ((fro->t <= (settime[f + 1] + GMX_REAL_EPS)) && do_enx(in, fr))
570 startstep_file = fr->step;
571 tadjust = settime[f] - fr->t;
572 if (cont_type[f + 1] == TIME_LAST)
574 settime[f + 1] = readtime[f + 1] - readtime[f] + settime[f];
575 cont_type[f + 1] = TIME_EXPLICIT;
580 if (tadjust + fr->t <= last_t)
582 /* Skip this frame, since we already have it / past it */
585 fprintf(debug, "fr->step %s, fr->t %.4f\n", gmx_step_str(fr->step, buf), fr->t);
586 fprintf(debug, "tadjust %12.6e + fr->t %12.6e <= t %12.6e\n", tadjust, fr->t, last_t);
591 fro->step = lastfilestep + fr->step - startstep_file;
592 fro->t = tadjust + fr->t;
594 bWrite = ((begin < 0 || (fro->t >= begin - GMX_REAL_EPS))
595 && (end < 0 || (fro->t <= end + GMX_REAL_EPS))
596 && (fro->t <= settime[f + 1] + 0.5 * timestep));
600 fprintf(debug, "fr->step %s, fr->t %.4f, fro->step %s fro->t %.4f, w %s\n",
601 gmx_step_str(fr->step, buf), fr->t, gmx_step_str(fro->step, buf2), fro->t,
602 gmx::boolToString(bWrite));
607 if ((end > 0) && (fro->t > end + GMX_REAL_EPS))
614 if (fro->t >= begin - GMX_REAL_EPS)
622 update_ee_sum(nre, &ee_sum_step, &ee_sum_nsteps, &ee_sum_nsum, ee_sum, fr, fro->step);
626 /* determine if we should write it */
627 if (bWrite && (delta_t == 0 || bRmod(fro->t, toffset, delta_t)))
629 laststep = fro->step;
634 fprintf(stderr, "\nContinue writing frames from t=%g, step=%s\n", fro->t,
635 gmx_step_str(fro->step, buf));
638 /* Copy the energies */
639 for (i = 0; i < nre; i++)
641 fro->ener[i].e = fr->ener[i].e;
644 fro->nsteps = ee_sum_nsteps;
647 if (ee_sum_nsum <= 1)
653 fro->nsum = int64_to_int(ee_sum_nsum, "energy average summation");
654 /* Copy the energy sums */
655 for (i = 0; i < nre; i++)
657 fro->ener[i].esum = ee_sum[i].esum;
658 fro->ener[i].eav = ee_sum[i].eav;
661 /* We wrote the energies, so reset the counts */
667 for (kkk = 0; kkk < nset; kkk++)
669 fro->ener[set[kkk]].e *= scalefac;
672 fro->ener[set[kkk]].eav *= scalefac * scalefac;
673 fro->ener[set[kkk]].esum *= scalefac;
677 /* Copy restraint stuff */
678 /*fro->ndisre = fr->ndisre;
679 fro->disre_rm3tav = fr->disre_rm3tav;
680 fro->disre_rt = fr->disre_rt;*/
681 fro->nblock = fr->nblock;
682 /*fro->nr = fr->nr;*/
683 fro->block = fr->block;
685 /* check if we have blocks with delta_h data and are throwing
692 if (!blocks || nblocks_alloc < fr->nblock)
694 /* we pre-allocate the blocks */
695 nblocks_alloc = fr->nblock;
696 snew(blocks, nblocks_alloc);
698 nblocks = 0; /* number of blocks so far */
700 for (i = 0; i < fr->nblock; i++)
702 if ((fr->block[i].id != enxDHCOLL) && (fr->block[i].id != enxDH)
703 && (fr->block[i].id != enxDHHIST))
705 /* copy everything verbatim */
706 blocks[nblocks] = fr->block[i];
710 /* now set the block pointer to the new blocks */
711 fro->nblock = nblocks;
714 else if (delta_t > 0)
716 if (!warned_about_dh)
718 for (i = 0; i < fr->nblock; i++)
720 if (fr->block[i].id == enxDH || fr->block[i].id == enxDHHIST)
723 if (fr->block[i].id == enxDH)
725 size = fr->block[i].sub[2].nr;
733 printf("\nWARNING: %s contains delta H blocks or "
734 "histograms for which\n"
735 " some data is thrown away on a "
736 "block-by-block basis, where each block\n"
737 " contains up to %d samples.\n"
738 " This is almost certainly not what you "
740 " Use the -rmdh option to throw all delta H "
742 " Use g_energy -odh option to extract these "
744 files[f].c_str(), size);
745 warned_about_dh = TRUE;
755 if (noutfr % 1000 == 0)
757 fprintf(stderr, "Writing frame time %g ", fro->t);
762 if (f == files.size())
766 printf("\nLast step written from %s: t %g, step %s\n", files[f].c_str(), last_t,
767 gmx_step_str(laststep, buf));
768 lastfilestep = laststep;
770 /* set the next time from the last in previous file */
771 if (cont_type[f + 1] == TIME_CONTINUE)
773 settime[f + 1] = fro->t;
774 /* in this case we have already written the last frame of
775 * previous file, so update begin to avoid doubling it
776 * with the start of the next file
778 begin = fro->t + 0.5 * timestep;
779 /* cont_type[f+1]==TIME_EXPLICIT; */
782 if ((fro->t < end) && (f < files.size() - 1) && (fro->t < settime[f + 1] - 1.5 * timestep))
784 fprintf(stderr, "\nWARNING: There might be a gap around t=%g\n", fro->t);
787 /* move energies to lastee */
789 free_enxnms(this_nre, enm);
791 fprintf(stderr, "\n");
795 fprintf(stderr, "No frames written.\n");
799 fprintf(stderr, "Last frame written was at step %s, time %f\n",
800 gmx_step_str(fro->step, buf), fro->t);
801 fprintf(stderr, "Wrote %d frames\n", noutfr);