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, 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.
43 #include "gromacs/legacyheaders/typedefs.h"
44 #include "gromacs/utility/smalloc.h"
45 #include "gromacs/commandline/pargs.h"
46 #include "gromacs/legacyheaders/disre.h"
47 #include "gromacs/legacyheaders/names.h"
48 #include "gromacs/legacyheaders/macros.h"
49 #include "gromacs/utility/fatalerror.h"
50 #include "gromacs/fileio/enxio.h"
51 #include "gromacs/math/vec.h"
53 #include "gromacs/fileio/trxio.h"
55 #define TIME_EXPLICIT 0
56 #define TIME_CONTINUE 1
62 static int *select_it(int nre, gmx_enxnm_t *nm, int *nset)
67 gmx_bool bVerbose = TRUE;
69 if ((getenv("GMX_ENER_VERBOSE")) != NULL)
74 fprintf(stderr, "Select the terms you want to scale from the following list\n");
75 fprintf(stderr, "End your selection with 0\n");
79 for (k = 0; (k < nre); )
81 for (j = 0; (j < 4) && (k < nre); j++, k++)
83 fprintf(stderr, " %3d=%14s", k+1, nm[k].name);
85 fprintf(stderr, "\n");
92 if (1 != scanf("%d", &n))
94 gmx_fatal(FARGS, "Cannot read energy term");
96 if ((n > 0) && (n <= nre))
104 for (i = (*nset) = 0; (i < nre); i++)
117 static void sort_files(char **fnms, real *settime, int nfile)
123 for (i = 0; i < nfile; i++)
126 for (j = i+1; j < nfile; j++)
128 if (settime[j] < settime[minidx])
135 timeswap = settime[i];
136 settime[i] = settime[minidx];
137 settime[minidx] = timeswap;
139 fnms[i] = fnms[minidx];
140 fnms[minidx] = chptr;
146 static int scan_ene_files(char **fnms, int nfiles,
147 real *readtime, real *timestep, int *nremax)
149 /* Check number of energy terms and start time of all files */
150 int f, i, nre, nremin = 0, nresav = 0;
153 char inputstring[STRLEN];
159 for (f = 0; f < nfiles; f++)
161 in = open_enx(fnms[f], "r");
163 do_enxnms(in, &nre, &enm);
180 nremin = min(nremin, fr->nre);
181 *nremax = max(*nremax, fr->nre);
185 "Energy files don't match, different number of energies:\n"
186 " %s: %d\n %s: %d\n", fnms[f-1], nresav, fnms[f], fr->nre);
188 "\nContinue conversion using only the first %d terms (n/y)?\n"
189 "(you should be sure that the energy terms match)\n", nremin);
190 if (NULL == fgets(inputstring, STRLEN-1, stdin))
192 gmx_fatal(FARGS, "Error reading user input");
194 if (inputstring[0] != 'y' && inputstring[0] != 'Y')
196 fprintf(stderr, "Will not convert\n");
205 fprintf(stderr, "\n");
206 free_enxnms(nre, enm);
216 static void edit_files(char **fnms, int nfiles, real *readtime,
217 real *settime, int *cont_type, gmx_bool bSetTime, gmx_bool bSort)
221 char inputstring[STRLEN], *chptr;
227 fprintf(stderr, "\n\nEnter the new start time:\n\n");
231 fprintf(stderr, "\n\nEnter the new start time for each file.\n"
232 "There are two special options, both disables sorting:\n\n"
233 "c (continue) - The start time is taken from the end\n"
234 "of the previous file. Use it when your continuation run\n"
235 "restarts with t=0 and there is no overlap.\n\n"
236 "l (last) - The time in this file will be changed the\n"
237 "same amount as in the previous. Use it when the time in the\n"
238 "new run continues from the end of the previous one,\n"
239 "since this takes possible overlap into account.\n\n");
242 fprintf(stderr, " File Current start New start\n"
243 "---------------------------------------------------------\n");
245 for (i = 0; i < nfiles; i++)
247 fprintf(stderr, "%25s %10.3f ", fnms[i], readtime[i]);
251 if (NULL == fgets(inputstring, STRLEN-1, stdin))
253 gmx_fatal(FARGS, "Error reading user input");
255 inputstring[strlen(inputstring)-1] = 0;
257 if (inputstring[0] == 'c' || inputstring[0] == 'C')
259 cont_type[i] = TIME_CONTINUE;
262 settime[i] = FLT_MAX;
264 else if (inputstring[0] == 'l' ||
265 inputstring[0] == 'L')
267 cont_type[i] = TIME_LAST;
270 settime[i] = FLT_MAX;
274 settime[i] = strtod(inputstring, &chptr);
275 if (chptr == inputstring)
277 fprintf(stderr, "Try that again: ");
281 cont_type[i] = TIME_EXPLICIT;
288 if (cont_type[0] != TIME_EXPLICIT)
290 cont_type[0] = TIME_EXPLICIT;
296 for (i = 0; i < nfiles; i++)
298 settime[i] = readtime[i];
302 if (bSort && (nfiles > 1))
304 sort_files(fnms, settime, nfiles);
308 fprintf(stderr, "Sorting disabled.\n");
312 /* Write out the new order and start times */
313 fprintf(stderr, "\nSummary of files and start times used:\n\n"
315 "-----------------------------------------\n");
316 for (i = 0; i < nfiles; i++)
318 switch (cont_type[i])
321 fprintf(stderr, "%25s %10.3f\n", fnms[i], settime[i]);
324 fprintf(stderr, "%25s Continue from end of last file\n", fnms[i]);
327 fprintf(stderr, "%25s Change by same amount as last file\n", fnms[i]);
331 fprintf(stderr, "\n");
333 settime[nfiles] = FLT_MAX;
334 cont_type[nfiles] = TIME_EXPLICIT;
335 readtime[nfiles] = FLT_MAX;
339 static void copy_ee(t_energy *src, t_energy *dst, int nre)
343 for (i = 0; i < nre; i++)
346 dst[i].esum = src[i].esum;
347 dst[i].eav = src[i].eav;
351 static void update_ee(t_energy *lastee, gmx_int64_t laststep,
352 t_energy *startee, gmx_int64_t startstep,
353 t_energy *ee, int step,
354 t_energy *outee, int nre)
357 double sigmacorr, nom, denom;
358 double prestart_esum;
359 double prestart_sigma;
361 for (i = 0; i < nre; i++)
363 outee[i].e = ee[i].e;
364 /* add statistics from earlier file if present */
367 outee[i].esum = lastee[i].esum+ee[i].esum;
368 nom = (lastee[i].esum*(step+1)-ee[i].esum*(laststep));
369 denom = ((step+1.0)*(laststep)*(step+1.0+laststep));
370 sigmacorr = nom*nom/denom;
371 outee[i].eav = lastee[i].eav+ee[i].eav+sigmacorr;
375 /* otherwise just copy to output */
376 outee[i].esum = ee[i].esum;
377 outee[i].eav = ee[i].eav;
380 /* if we didnt start to write at the first frame
381 * we must compensate the statistics for this
382 * there are laststep frames in the earlier file
383 * and step+1 frames in this one.
387 gmx_int64_t q = laststep+step;
388 gmx_int64_t p = startstep+1;
389 prestart_esum = startee[i].esum-startee[i].e;
390 sigmacorr = prestart_esum-(p-1)*startee[i].e;
391 prestart_sigma = startee[i].eav-
392 sigmacorr*sigmacorr/(p*(p-1));
393 sigmacorr = prestart_esum/(p-1)-
395 outee[i].esum -= prestart_esum;
398 outee[i].eav = outee[i].eav-prestart_sigma-
399 sigmacorr*sigmacorr*((p-1)*q)/(q-p+1);
403 if ((outee[i].eav/(laststep+step+1)) < (GMX_REAL_EPS))
410 static void update_ee_sum(int nre,
411 gmx_int64_t *ee_sum_step,
412 gmx_int64_t *ee_sum_nsteps,
413 gmx_int64_t *ee_sum_nsum,
415 t_enxframe *fr, int out_step)
417 gmx_int64_t nsteps, nsum, fr_nsum;
420 nsteps = *ee_sum_nsteps;
433 for (i = 0; i < nre; i++)
435 ee_sum[i].esum = fr->ener[i].e;
441 for (i = 0; i < nre; i++)
443 ee_sum[i].esum = fr->ener[i].esum;
444 ee_sum[i].eav = fr->ener[i].eav;
450 else if (out_step + *ee_sum_nsum - *ee_sum_step == nsteps + fr->nsteps)
454 for (i = 0; i < nre; i++)
457 dsqr(ee_sum[i].esum/nsum
458 - (ee_sum[i].esum + fr->ener[i].e)/(nsum + 1))*nsum*(nsum + 1);
459 ee_sum[i].esum += fr->ener[i].e;
464 for (i = 0; i < fr->nre; i++)
468 dsqr(ee_sum[i].esum/nsum
469 - (ee_sum[i].esum + fr->ener[i].esum)/(nsum + fr->nsum))*
470 nsum*(nsum + fr->nsum)/(double)fr->nsum;
471 ee_sum[i].esum += fr->ener[i].esum;
474 nsteps += fr->nsteps;
481 fprintf(stderr, "\nWARNING: missing energy sums at time %f\n", fr->t);
487 *ee_sum_step = out_step;
488 *ee_sum_nsteps = nsteps;
492 int gmx_eneconv(int argc, char *argv[])
494 const char *desc[] = {
495 "With [IT]multiple files[it] specified for the [TT]-f[tt] option:[BR]",
496 "Concatenates several energy files in sorted order.",
497 "In the case of double time frames, the one",
498 "in the later file is used. By specifying [TT]-settime[tt] you will be",
499 "asked for the start time of each file. The input files are taken",
500 "from the command line,",
501 "such that the command [TT]gmx eneconv -f *.edr -o fixed.edr[tt] should do",
503 "With [IT]one file[it] specified for [TT]-f[tt]:[BR]",
504 "Reads one energy file and writes another, applying the [TT]-dt[tt],",
505 "[TT]-offset[tt], [TT]-t0[tt] and [TT]-settime[tt] options and",
506 "converting to a different format if necessary (indicated by file",
508 "[TT]-settime[tt] is applied first, then [TT]-dt[tt]/[TT]-offset[tt]",
509 "followed by [TT]-b[tt] and [TT]-e[tt] to select which frames to write."
511 const char *bugs[] = {
512 "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."
514 ener_file_t in = NULL, out = NULL;
515 gmx_enxnm_t *enm = NULL;
517 ener_file_t in, out = NULL;
518 gmx_enxnm_t *enm = NULL;
520 t_enxframe *fr, *fro;
521 gmx_int64_t ee_sum_step = 0, ee_sum_nsteps, ee_sum_nsum;
523 gmx_int64_t lastfilestep, laststep, startstep, startstep_file = 0;
525 int nre, nremax, this_nre, nfile, f, i, j, kkk, nset, *set = NULL;
528 real *readtime, *settime, timestep, t1, tadjust;
529 char inputstring[STRLEN], *chptr, buf[22], buf2[22], buf3[22];
532 gmx_bool bNewFile, bFirst, bNewOutput;
534 gmx_bool warned_about_dh = FALSE;
535 t_enxblock *blocks = NULL;
537 int nblocks_alloc = 0;
540 { efEDR, "-f", NULL, ffRDMULT },
541 { efEDR, "-o", "fixed", ffWRITE },
544 #define NFILE asize(fnm)
546 static real delta_t = 0.0, toffset = 0, scalefac = 1;
547 static gmx_bool bSetTime = FALSE;
548 static gmx_bool bSort = TRUE, bError = TRUE;
549 static real begin = -1;
550 static real end = -1;
551 gmx_bool remove_dh = FALSE;
554 { "-b", FALSE, etREAL, {&begin},
555 "First time to use"},
556 { "-e", FALSE, etREAL, {&end},
558 { "-dt", FALSE, etREAL, {&delta_t},
559 "Only write out frame when t MOD dt = offset" },
560 { "-offset", FALSE, etREAL, {&toffset},
561 "Time offset for [TT]-dt[tt] option" },
562 { "-settime", FALSE, etBOOL, {&bSetTime},
563 "Change starting time interactively" },
564 { "-sort", FALSE, etBOOL, {&bSort},
565 "Sort energy files (not frames)"},
566 { "-rmdh", FALSE, etBOOL, {&remove_dh},
567 "Remove free energy block data" },
568 { "-scalefac", FALSE, etREAL, {&scalefac},
569 "Multiply energy component by this factor" },
570 { "-error", FALSE, etBOOL, {&bError},
571 "Stop on errors in the file" }
574 if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa),
575 pa, asize(desc), desc, asize(bugs), bugs, &oenv))
585 laststep = startstep = 0;
587 nfile = opt2fns(&fnms, "-f", NFILE, fnm);
591 gmx_fatal(FARGS, "No input files!");
594 snew(settime, nfile+1);
595 snew(readtime, nfile+1);
596 snew(cont_type, nfile+1);
598 nre = scan_ene_files(fnms, nfile, readtime, ×tep, &nremax);
599 edit_files(fnms, nfile, readtime, settime, cont_type, bSetTime, bSort);
603 snew(ee_sum, nremax);
609 snew(fro->ener, nremax);
615 for (f = 0; f < nfile; f++)
619 in = open_enx(fnms[f], "r");
621 do_enxnms(in, &this_nre, &enm);
626 set = select_it(nre, enm, &nset);
629 /* write names to the output file */
630 out = open_enx(opt2fn("-o", NFILE, fnm), "w");
631 do_enxnms(out, &nre, &enm);
634 /* start reading from the next file */
635 while ((fro->t <= (settime[f+1] + GMX_REAL_EPS)) &&
640 startstep_file = fr->step;
641 tadjust = settime[f] - fr->t;
642 if (cont_type[f+1] == TIME_LAST)
644 settime[f+1] = readtime[f+1]-readtime[f]+settime[f];
645 cont_type[f+1] = TIME_EXPLICIT;
650 if (tadjust + fr->t <= last_t)
652 /* Skip this frame, since we already have it / past it */
655 fprintf(debug, "fr->step %s, fr->t %.4f\n",
656 gmx_step_str(fr->step, buf), fr->t);
657 fprintf(debug, "tadjust %12.6e + fr->t %12.6e <= t %12.6e\n",
658 tadjust, fr->t, last_t);
663 fro->step = lastfilestep + fr->step - startstep_file;
664 fro->t = tadjust + fr->t;
666 bWrite = ((begin < 0 || (begin >= 0 && (fro->t >= begin-GMX_REAL_EPS))) &&
667 (end < 0 || (end >= 0 && (fro->t <= end +GMX_REAL_EPS))) &&
668 (fro->t <= settime[f+1]+0.5*timestep));
673 "fr->step %s, fr->t %.4f, fro->step %s fro->t %.4f, w %d\n",
674 gmx_step_str(fr->step, buf), fr->t,
675 gmx_step_str(fro->step, buf2), fro->t, bWrite);
680 if ((end > 0) && (fro->t > end+GMX_REAL_EPS))
687 if (fro->t >= begin-GMX_REAL_EPS)
692 startstep = fr->step;
696 update_ee_sum(nre, &ee_sum_step, &ee_sum_nsteps, &ee_sum_nsum, ee_sum,
701 /* determine if we should write it */
702 if (bWrite && (delta_t == 0 || bRmod(fro->t, toffset, delta_t)))
704 laststep = fro->step;
709 fprintf(stderr, "\nContinue writing frames from t=%g, step=%s\n",
710 fro->t, gmx_step_str(fro->step, buf));
713 /* Copy the energies */
714 for (i = 0; i < nre; i++)
716 fro->ener[i].e = fr->ener[i].e;
719 fro->nsteps = ee_sum_nsteps;
722 if (ee_sum_nsum <= 1)
728 fro->nsum = gmx_int64_to_int(ee_sum_nsum,
729 "energy average summation");
730 /* Copy the energy sums */
731 for (i = 0; i < nre; i++)
733 fro->ener[i].esum = ee_sum[i].esum;
734 fro->ener[i].eav = ee_sum[i].eav;
737 /* We wrote the energies, so reset the counts */
743 for (kkk = 0; kkk < nset; kkk++)
745 fro->ener[set[kkk]].e *= scalefac;
748 fro->ener[set[kkk]].eav *= scalefac*scalefac;
749 fro->ener[set[kkk]].esum *= scalefac;
753 /* Copy restraint stuff */
754 /*fro->ndisre = fr->ndisre;
755 fro->disre_rm3tav = fr->disre_rm3tav;
756 fro->disre_rt = fr->disre_rt;*/
757 fro->nblock = fr->nblock;
758 /*fro->nr = fr->nr;*/
759 fro->block = fr->block;
761 /* check if we have blocks with delta_h data and are throwing
768 if (!blocks || nblocks_alloc < fr->nblock)
770 /* we pre-allocate the blocks */
771 nblocks_alloc = fr->nblock;
772 snew(blocks, nblocks_alloc);
774 nblocks = 0; /* number of blocks so far */
776 for (i = 0; i < fr->nblock; i++)
778 if ( (fr->block[i].id != enxDHCOLL) &&
779 (fr->block[i].id != enxDH) &&
780 (fr->block[i].id != enxDHHIST) )
782 /* copy everything verbatim */
783 blocks[nblocks] = fr->block[i];
787 /* now set the block pointer to the new blocks */
788 fro->nblock = nblocks;
791 else if (delta_t > 0)
793 if (!warned_about_dh)
795 for (i = 0; i < fr->nblock; i++)
797 if (fr->block[i].id == enxDH ||
798 fr->block[i].id == enxDHHIST)
801 if (fr->block[i].id == enxDH)
803 size = fr->block[i].sub[2].nr;
811 printf("\nWARNING: %s contains delta H blocks or histograms for which\n"
812 " some data is thrown away on a block-by-block basis, where each block\n"
813 " contains up to %d samples.\n"
814 " This is almost certainly not what you want.\n"
815 " Use the -rmdh option to throw all delta H samples away.\n"
816 " Use g_energy -odh option to extract these samples.\n",
818 warned_about_dh = TRUE;
828 if (noutfr % 1000 == 0)
830 fprintf(stderr, "Writing frame time %g ", fro->t);
839 printf("\nLast step written from %s: t %g, step %s\n",
840 fnms[f], last_t, gmx_step_str(laststep, buf));
841 lastfilestep = laststep;
843 /* set the next time from the last in previous file */
844 if (cont_type[f+1] == TIME_CONTINUE)
846 settime[f+1] = fro->t;
847 /* in this case we have already written the last frame of
848 * previous file, so update begin to avoid doubling it
849 * with the start of the next file
851 begin = fro->t+0.5*timestep;
852 /* cont_type[f+1]==TIME_EXPLICIT; */
855 if ((fro->t < end) && (f < nfile-1) &&
856 (fro->t < settime[f+1]-1.5*timestep))
859 "\nWARNING: There might be a gap around t=%g\n", fro->t);
862 /* move energies to lastee */
864 free_enxnms(this_nre, enm);
866 fprintf(stderr, "\n");
870 fprintf(stderr, "No frames written.\n");
874 fprintf(stderr, "Last frame written was at step %s, time %f\n",
875 gmx_step_str(fro->step, buf), fro->t);
876 fprintf(stderr, "Wrote %d frames\n", noutfr);