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.
45 #include "gromacs/commandline/pargs.h"
46 #include "gromacs/fileio/enxio.h"
47 #include "gromacs/fileio/trxio.h"
48 #include "gromacs/gmxana/gmx_ana.h"
49 #include "gromacs/listed_forces/disre.h"
50 #include "gromacs/math/functions.h"
51 #include "gromacs/math/vec.h"
52 #include "gromacs/mdtypes/md_enums.h"
53 #include "gromacs/trajectory/energyframe.h"
54 #include "gromacs/utility/arrayref.h"
55 #include "gromacs/utility/arraysize.h"
56 #include "gromacs/utility/cstringutil.h"
57 #include "gromacs/utility/fatalerror.h"
58 #include "gromacs/utility/int64_to_int.h"
59 #include "gromacs/utility/smalloc.h"
60 #include "gromacs/utility/strconvert.h"
62 #define TIME_EXPLICIT 0
63 #define TIME_CONTINUE 1
69 static int *select_it(int nre, gmx_enxnm_t *nm, int *nset)
74 gmx_bool bVerbose = TRUE;
76 if ((getenv("GMX_ENER_VERBOSE")) != nullptr)
81 fprintf(stderr, "Select the terms you want to scale from the following list\n");
82 fprintf(stderr, "End your selection with 0\n");
86 for (k = 0; (k < nre); )
88 for (j = 0; (j < 4) && (k < nre); j++, k++)
90 fprintf(stderr, " %3d=%14s", k+1, nm[k].name);
92 fprintf(stderr, "\n");
99 if (1 != scanf("%d", &n))
101 gmx_fatal(FARGS, "Cannot read energy term");
103 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,
150 real *readtime, real *timestep, int *nremax)
152 /* Check number of energy terms and start time of all files */
153 int nre, nremin = 0, nresav = 0;
156 char inputstring[STRLEN];
162 for (size_t f = 0; f < files.size(); f++)
164 in = open_enx(files[f].c_str(), "r");
166 do_enxnms(in, &nre, &enm);
183 nremin = std::min(nremin, fr->nre);
184 *nremax = std::max(*nremax, fr->nre);
188 "Energy files don't match, different number of energies:\n"
189 " %s: %d\n %s: %d\n", 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", nremin);
193 if (nullptr == fgets(inputstring, STRLEN-1, stdin))
195 gmx_fatal(FARGS, "Error reading user input");
197 if (inputstring[0] != 'y' && inputstring[0] != 'Y')
199 fprintf(stderr, "Will not convert\n");
208 fprintf(stderr, "\n");
209 free_enxnms(nre, enm);
219 static void edit_files(gmx::ArrayRef<std::string> files, real *readtime,
220 real *settime, int *cont_type, gmx_bool bSetTime, gmx_bool bSort)
223 char inputstring[STRLEN], *chptr;
227 if (files.size() == 1)
229 fprintf(stderr, "\n\nEnter the new start time:\n\n");
233 fprintf(stderr, "\n\nEnter the new start time for each file.\n"
234 "There are two special options, both disables sorting:\n\n"
235 "c (continue) - The start time is taken from the end\n"
236 "of the previous file. Use it when your continuation run\n"
237 "restarts with t=0 and there is no overlap.\n\n"
238 "l (last) - The time in this file will be changed the\n"
239 "same amount as in the previous. Use it when the time in the\n"
240 "new run continues from the end of the previous one,\n"
241 "since this takes possible overlap into account.\n\n");
244 fprintf(stderr, " File Current start New start\n"
245 "---------------------------------------------------------\n");
247 for (gmx::index i = 0; i < files.ssize(); i++)
249 fprintf(stderr, "%25s %10.3f ", files[i].c_str(), readtime[i]);
253 if (nullptr == fgets(inputstring, STRLEN-1, stdin))
255 gmx_fatal(FARGS, "Error reading user input");
257 inputstring[std::strlen(inputstring)-1] = 0;
259 if (inputstring[0] == 'c' || inputstring[0] == 'C')
261 cont_type[i] = TIME_CONTINUE;
264 settime[i] = FLT_MAX;
266 else if (inputstring[0] == 'l' ||
267 inputstring[0] == 'L')
269 cont_type[i] = TIME_LAST;
272 settime[i] = FLT_MAX;
276 settime[i] = strtod(inputstring, &chptr);
277 if (chptr == inputstring)
279 fprintf(stderr, "Try that again: ");
283 cont_type[i] = TIME_EXPLICIT;
290 if (cont_type[0] != TIME_EXPLICIT)
292 cont_type[0] = TIME_EXPLICIT;
298 for (gmx::index i = 0; i < files.ssize(); i++)
300 settime[i] = readtime[i];
304 if (bSort && files.size() > 1)
306 sort_files(files, settime);
310 fprintf(stderr, "Sorting disabled.\n");
314 /* Write out the new order and start times */
315 fprintf(stderr, "\nSummary of files and start times used:\n\n"
317 "-----------------------------------------\n");
318 for (gmx::index i = 0; i < files.ssize(); i++)
320 switch (cont_type[i])
323 fprintf(stderr, "%25s %10.3f\n", files[i].c_str(), settime[i]);
326 fprintf(stderr, "%25s Continue from end of last file\n", files[i].c_str());
329 fprintf(stderr, "%25s Change by same amount as last file\n", files[i].c_str());
333 fprintf(stderr, "\n");
335 settime[files.size()] = FLT_MAX;
336 cont_type[files.size()] = TIME_EXPLICIT;
337 readtime[files.size()] = FLT_MAX;
341 static void update_ee_sum(int nre,
342 int64_t *ee_sum_step,
343 int64_t *ee_sum_nsteps,
344 int64_t *ee_sum_nsum,
346 t_enxframe *fr, int out_step)
348 int64_t nsteps, nsum, fr_nsum;
351 nsteps = *ee_sum_nsteps;
364 for (i = 0; i < nre; i++)
366 ee_sum[i].esum = fr->ener[i].e;
372 for (i = 0; i < nre; i++)
374 ee_sum[i].esum = fr->ener[i].esum;
375 ee_sum[i].eav = fr->ener[i].eav;
381 else if (out_step + *ee_sum_nsum - *ee_sum_step == nsteps + fr->nsteps)
385 for (i = 0; i < nre; i++)
388 gmx::square(ee_sum[i].esum/nsum
389 - (ee_sum[i].esum + fr->ener[i].e)/(nsum + 1))*nsum*(nsum + 1);
390 ee_sum[i].esum += fr->ener[i].e;
395 for (i = 0; i < fr->nre; i++)
399 gmx::square(ee_sum[i].esum/nsum
400 - (ee_sum[i].esum + fr->ener[i].esum)/(nsum + fr->nsum))*
401 nsum*(nsum + fr->nsum)/static_cast<double>(fr->nsum);
402 ee_sum[i].esum += fr->ener[i].esum;
405 nsteps += fr->nsteps;
412 fprintf(stderr, "\nWARNING: missing energy sums at time %f\n", fr->t);
418 *ee_sum_step = out_step;
419 *ee_sum_nsteps = nsteps;
423 int gmx_eneconv(int argc, char *argv[])
425 const char *desc[] = {
426 "With [IT]multiple files[it] specified for the [TT]-f[tt] option:[PAR]",
427 "Concatenates several energy files in sorted order.",
428 "In the case of double time frames, the one",
429 "in the later file is used. By specifying [TT]-settime[tt] you will be",
430 "asked for the start time of each file. The input files are taken",
431 "from the command line,",
432 "such that the command [TT]gmx eneconv -f *.edr -o fixed.edr[tt] should do",
434 "With [IT]one file[it] specified for [TT]-f[tt]:[PAR]",
435 "Reads one energy file and writes another, applying the [TT]-dt[tt],",
436 "[TT]-offset[tt], [TT]-t0[tt] and [TT]-settime[tt] options and",
437 "converting to a different format if necessary (indicated by file",
439 "[TT]-settime[tt] is applied first, then [TT]-dt[tt]/[TT]-offset[tt]",
440 "followed by [TT]-b[tt] and [TT]-e[tt] to select which frames to write."
442 const char *bugs[] = {
443 "When combining trajectories the sigma and E^2 (necessary for statistics) are not updated correctly. Only the actual energy is correct. One thus has to compute statistics in another way."
445 ener_file_t in = nullptr, out = nullptr;
446 gmx_enxnm_t *enm = nullptr;
448 ener_file_t in, out = NULL;
449 gmx_enxnm_t *enm = NULL;
451 t_enxframe *fr, *fro;
452 int64_t ee_sum_step = 0, ee_sum_nsteps, ee_sum_nsum;
454 int64_t lastfilestep, laststep, startstep_file = 0;
456 int nre, nremax, this_nre, i, kkk, nset, *set = nullptr;
458 real *readtime, *settime, timestep, tadjust;
459 char buf[22], buf2[22];
461 gmx_bool bNewFile, bFirst, bNewOutput;
462 gmx_output_env_t *oenv;
463 gmx_bool warned_about_dh = FALSE;
464 t_enxblock *blocks = nullptr;
466 int nblocks_alloc = 0;
469 { efEDR, "-f", nullptr, ffRDMULT },
470 { efEDR, "-o", "fixed", ffWRITE },
473 #define NFILE asize(fnm)
475 static real delta_t = 0.0, toffset = 0, scalefac = 1;
476 static gmx_bool bSetTime = FALSE;
477 static gmx_bool bSort = TRUE, bError = TRUE;
478 static real begin = -1;
479 static real end = -1;
480 gmx_bool remove_dh = FALSE;
483 { "-b", FALSE, etREAL, {&begin},
484 "First time to use"},
485 { "-e", FALSE, etREAL, {&end},
487 { "-dt", FALSE, etREAL, {&delta_t},
488 "Only write out frame when t MOD dt = offset" },
489 { "-offset", FALSE, etREAL, {&toffset},
490 "Time offset for [TT]-dt[tt] option" },
491 { "-settime", FALSE, etBOOL, {&bSetTime},
492 "Change starting time interactively" },
493 { "-sort", FALSE, etBOOL, {&bSort},
494 "Sort energy files (not frames)"},
495 { "-rmdh", FALSE, etBOOL, {&remove_dh},
496 "Remove free energy block data" },
497 { "-scalefac", FALSE, etREAL, {&scalefac},
498 "Multiply energy component by this factor" },
499 { "-error", FALSE, etBOOL, {&bError},
500 "Stop on errors in the file" }
503 if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa),
504 pa, asize(desc), desc, asize(bugs), bugs, &oenv))
508 fprintf(stdout, "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)) &&
571 startstep_file = fr->step;
572 tadjust = settime[f] - fr->t;
573 if (cont_type[f+1] == TIME_LAST)
575 settime[f+1] = readtime[f+1]-readtime[f]+settime[f];
576 cont_type[f+1] = TIME_EXPLICIT;
581 if (tadjust + fr->t <= last_t)
583 /* Skip this frame, since we already have it / past it */
586 fprintf(debug, "fr->step %s, fr->t %.4f\n",
587 gmx_step_str(fr->step, buf), fr->t);
588 fprintf(debug, "tadjust %12.6e + fr->t %12.6e <= t %12.6e\n",
589 tadjust, fr->t, last_t);
594 fro->step = lastfilestep + fr->step - startstep_file;
595 fro->t = tadjust + fr->t;
597 bWrite = ((begin < 0 || (fro->t >= begin-GMX_REAL_EPS)) &&
598 (end < 0 || (fro->t <= end +GMX_REAL_EPS)) &&
599 (fro->t <= settime[f+1]+0.5*timestep));
604 "fr->step %s, fr->t %.4f, fro->step %s fro->t %.4f, w %s\n",
605 gmx_step_str(fr->step, buf), fr->t,
606 gmx_step_str(fro->step, buf2), fro->t, gmx::boolToString(bWrite));
611 if ((end > 0) && (fro->t > end+GMX_REAL_EPS))
618 if (fro->t >= begin-GMX_REAL_EPS)
626 update_ee_sum(nre, &ee_sum_step, &ee_sum_nsteps, &ee_sum_nsum, ee_sum,
631 /* determine if we should write it */
632 if (bWrite && (delta_t == 0 || bRmod(fro->t, toffset, delta_t)))
634 laststep = fro->step;
639 fprintf(stderr, "\nContinue writing frames from t=%g, step=%s\n",
640 fro->t, gmx_step_str(fro->step, buf));
643 /* Copy the energies */
644 for (i = 0; i < nre; i++)
646 fro->ener[i].e = fr->ener[i].e;
649 fro->nsteps = ee_sum_nsteps;
652 if (ee_sum_nsum <= 1)
658 fro->nsum = int64_to_int(ee_sum_nsum,
659 "energy average summation");
660 /* Copy the energy sums */
661 for (i = 0; i < nre; i++)
663 fro->ener[i].esum = ee_sum[i].esum;
664 fro->ener[i].eav = ee_sum[i].eav;
667 /* We wrote the energies, so reset the counts */
673 for (kkk = 0; kkk < nset; kkk++)
675 fro->ener[set[kkk]].e *= scalefac;
678 fro->ener[set[kkk]].eav *= scalefac*scalefac;
679 fro->ener[set[kkk]].esum *= scalefac;
683 /* Copy restraint stuff */
684 /*fro->ndisre = fr->ndisre;
685 fro->disre_rm3tav = fr->disre_rm3tav;
686 fro->disre_rt = fr->disre_rt;*/
687 fro->nblock = fr->nblock;
688 /*fro->nr = fr->nr;*/
689 fro->block = fr->block;
691 /* check if we have blocks with delta_h data and are throwing
698 if (!blocks || nblocks_alloc < fr->nblock)
700 /* we pre-allocate the blocks */
701 nblocks_alloc = fr->nblock;
702 snew(blocks, nblocks_alloc);
704 nblocks = 0; /* number of blocks so far */
706 for (i = 0; i < fr->nblock; i++)
708 if ( (fr->block[i].id != enxDHCOLL) &&
709 (fr->block[i].id != enxDH) &&
710 (fr->block[i].id != enxDHHIST) )
712 /* copy everything verbatim */
713 blocks[nblocks] = fr->block[i];
717 /* now set the block pointer to the new blocks */
718 fro->nblock = nblocks;
721 else if (delta_t > 0)
723 if (!warned_about_dh)
725 for (i = 0; i < fr->nblock; i++)
727 if (fr->block[i].id == enxDH ||
728 fr->block[i].id == enxDHHIST)
731 if (fr->block[i].id == enxDH)
733 size = fr->block[i].sub[2].nr;
741 printf("\nWARNING: %s contains delta H blocks or histograms for which\n"
742 " some data is thrown away on a block-by-block basis, where each block\n"
743 " contains up to %d samples.\n"
744 " This is almost certainly not what you want.\n"
745 " Use the -rmdh option to throw all delta H samples away.\n"
746 " Use g_energy -odh option to extract these samples.\n",
747 files[f].c_str(), size);
748 warned_about_dh = TRUE;
758 if (noutfr % 1000 == 0)
760 fprintf(stderr, "Writing frame time %g ", fro->t);
765 if (f == files.size())
769 printf("\nLast step written from %s: t %g, step %s\n",
770 files[f].c_str(), last_t, gmx_step_str(laststep, buf));
771 lastfilestep = laststep;
773 /* set the next time from the last in previous file */
774 if (cont_type[f+1] == TIME_CONTINUE)
776 settime[f+1] = fro->t;
777 /* in this case we have already written the last frame of
778 * previous file, so update begin to avoid doubling it
779 * with the start of the next file
781 begin = fro->t+0.5*timestep;
782 /* cont_type[f+1]==TIME_EXPLICIT; */
785 if ((fro->t < end) && (f < files.size() - 1) &&
786 (fro->t < settime[f+1]-1.5*timestep))
789 "\nWARNING: There might be a gap around t=%g\n", fro->t);
792 /* move energies to lastee */
794 free_enxnms(this_nre, enm);
796 fprintf(stderr, "\n");
800 fprintf(stderr, "No frames written.\n");
804 fprintf(stderr, "Last frame written was at step %s, time %f\n",
805 gmx_step_str(fro->step, buf), fro->t);
806 fprintf(stderr, "Wrote %d frames\n", noutfr);