3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * Green Red Orange Magenta Azure Cyan Skyblue
48 #include "gmx_fatal.h"
53 #define TIME_EXPLICIT 0
54 #define TIME_CONTINUE 1
60 static int *select_it(int nre, gmx_enxnm_t *nm, int *nset)
65 gmx_bool bVerbose = TRUE;
67 if ((getenv("VERBOSE")) != NULL)
72 fprintf(stderr, "Select the terms you want to scale from the following list\n");
73 fprintf(stderr, "End your selection with 0\n");
77 for (k = 0; (k < nre); )
79 for (j = 0; (j < 4) && (k < nre); j++, k++)
81 fprintf(stderr, " %3d=%14s", k+1, nm[k].name);
83 fprintf(stderr, "\n");
90 if (1 != scanf("%d", &n))
92 gmx_fatal(FARGS, "Cannot read energy term");
94 if ((n > 0) && (n <= nre))
102 for (i = (*nset) = 0; (i < nre); i++)
115 static void sort_files(char **fnms, real *settime, int nfile)
121 for (i = 0; i < nfile; i++)
124 for (j = i+1; j < nfile; j++)
126 if (settime[j] < settime[minidx])
133 timeswap = settime[i];
134 settime[i] = settime[minidx];
135 settime[minidx] = timeswap;
137 fnms[i] = fnms[minidx];
138 fnms[minidx] = chptr;
144 static int scan_ene_files(char **fnms, int nfiles,
145 real *readtime, real *timestep, int *nremax)
147 /* Check number of energy terms and start time of all files */
148 int f, i, nre, nremin = 0, nresav = 0;
151 char inputstring[STRLEN];
157 for (f = 0; f < nfiles; f++)
159 in = open_enx(fnms[f], "r");
161 do_enxnms(in, &nre, &enm);
178 nremin = min(nremin, fr->nre);
179 *nremax = max(*nremax, fr->nre);
183 "Energy files don't match, different number of energies:\n"
184 " %s: %d\n %s: %d\n", fnms[f-1], nresav, fnms[f], fr->nre);
186 "\nContinue conversion using only the first %d terms (n/y)?\n"
187 "(you should be sure that the energy terms match)\n", nremin);
188 if (NULL == fgets(inputstring, STRLEN-1, stdin))
190 gmx_fatal(FARGS, "Error reading user input");
192 if (inputstring[0] != 'y' && inputstring[0] != 'Y')
194 fprintf(stderr, "Will not convert\n");
203 fprintf(stderr, "\n");
204 free_enxnms(nre, enm);
214 static void edit_files(char **fnms, int nfiles, real *readtime,
215 real *settime, int *cont_type, gmx_bool bSetTime, gmx_bool bSort)
219 char inputstring[STRLEN], *chptr;
225 fprintf(stderr, "\n\nEnter the new start time:\n\n");
229 fprintf(stderr, "\n\nEnter the new start time for each file.\n"
230 "There are two special options, both disables sorting:\n\n"
231 "c (continue) - The start time is taken from the end\n"
232 "of the previous file. Use it when your continuation run\n"
233 "restarts with t=0 and there is no overlap.\n\n"
234 "l (last) - The time in this file will be changed the\n"
235 "same amount as in the previous. Use it when the time in the\n"
236 "new run continues from the end of the previous one,\n"
237 "since this takes possible overlap into account.\n\n");
240 fprintf(stderr, " File Current start New start\n"
241 "---------------------------------------------------------\n");
243 for (i = 0; i < nfiles; i++)
245 fprintf(stderr, "%25s %10.3f ", fnms[i], readtime[i]);
249 if (NULL == fgets(inputstring, STRLEN-1, stdin))
251 gmx_fatal(FARGS, "Error reading user input");
253 inputstring[strlen(inputstring)-1] = 0;
255 if (inputstring[0] == 'c' || inputstring[0] == 'C')
257 cont_type[i] = TIME_CONTINUE;
260 settime[i] = FLT_MAX;
262 else if (inputstring[0] == 'l' ||
263 inputstring[0] == 'L')
265 cont_type[i] = TIME_LAST;
268 settime[i] = FLT_MAX;
272 settime[i] = strtod(inputstring, &chptr);
273 if (chptr == inputstring)
275 fprintf(stderr, "Try that again: ");
279 cont_type[i] = TIME_EXPLICIT;
286 if (cont_type[0] != TIME_EXPLICIT)
288 cont_type[0] = TIME_EXPLICIT;
294 for (i = 0; i < nfiles; i++)
296 settime[i] = readtime[i];
300 if (bSort && (nfiles > 1))
302 sort_files(fnms, settime, nfiles);
306 fprintf(stderr, "Sorting disabled.\n");
310 /* Write out the new order and start times */
311 fprintf(stderr, "\nSummary of files and start times used:\n\n"
313 "-----------------------------------------\n");
314 for (i = 0; i < nfiles; i++)
316 switch (cont_type[i])
319 fprintf(stderr, "%25s %10.3f\n", fnms[i], settime[i]);
322 fprintf(stderr, "%25s Continue from end of last file\n", fnms[i]);
325 fprintf(stderr, "%25s Change by same amount as last file\n", fnms[i]);
329 fprintf(stderr, "\n");
331 settime[nfiles] = FLT_MAX;
332 cont_type[nfiles] = TIME_EXPLICIT;
333 readtime[nfiles] = FLT_MAX;
337 static void copy_ee(t_energy *src, t_energy *dst, int nre)
341 for (i = 0; i < nre; i++)
344 dst[i].esum = src[i].esum;
345 dst[i].eav = src[i].eav;
349 static void update_ee(t_energy *lastee, gmx_large_int_t laststep,
350 t_energy *startee, gmx_large_int_t startstep,
351 t_energy *ee, int step,
352 t_energy *outee, int nre)
355 double sigmacorr, nom, denom;
356 double prestart_esum;
357 double prestart_sigma;
359 for (i = 0; i < nre; i++)
361 outee[i].e = ee[i].e;
362 /* add statistics from earlier file if present */
365 outee[i].esum = lastee[i].esum+ee[i].esum;
366 nom = (lastee[i].esum*(step+1)-ee[i].esum*(laststep));
367 denom = ((step+1.0)*(laststep)*(step+1.0+laststep));
368 sigmacorr = nom*nom/denom;
369 outee[i].eav = lastee[i].eav+ee[i].eav+sigmacorr;
373 /* otherwise just copy to output */
374 outee[i].esum = ee[i].esum;
375 outee[i].eav = ee[i].eav;
378 /* if we didnt start to write at the first frame
379 * we must compensate the statistics for this
380 * there are laststep frames in the earlier file
381 * and step+1 frames in this one.
385 gmx_large_int_t q = laststep+step;
386 gmx_large_int_t p = startstep+1;
387 prestart_esum = startee[i].esum-startee[i].e;
388 sigmacorr = prestart_esum-(p-1)*startee[i].e;
389 prestart_sigma = startee[i].eav-
390 sigmacorr*sigmacorr/(p*(p-1));
391 sigmacorr = prestart_esum/(p-1)-
393 outee[i].esum -= prestart_esum;
396 outee[i].eav = outee[i].eav-prestart_sigma-
397 sigmacorr*sigmacorr*((p-1)*q)/(q-p+1);
401 if ((outee[i].eav/(laststep+step+1)) < (GMX_REAL_EPS))
408 static void update_ee_sum(int nre,
409 gmx_large_int_t *ee_sum_step,
410 gmx_large_int_t *ee_sum_nsteps,
411 gmx_large_int_t *ee_sum_nsum,
413 t_enxframe *fr, int out_step)
415 gmx_large_int_t nsteps, nsum, fr_nsum;
418 nsteps = *ee_sum_nsteps;
431 for (i = 0; i < nre; i++)
433 ee_sum[i].esum = fr->ener[i].e;
439 for (i = 0; i < nre; i++)
441 ee_sum[i].esum = fr->ener[i].esum;
442 ee_sum[i].eav = fr->ener[i].eav;
448 else if (out_step + *ee_sum_nsum - *ee_sum_step == nsteps + fr->nsteps)
452 for (i = 0; i < nre; i++)
455 dsqr(ee_sum[i].esum/nsum
456 - (ee_sum[i].esum + fr->ener[i].e)/(nsum + 1))*nsum*(nsum + 1);
457 ee_sum[i].esum += fr->ener[i].e;
462 for (i = 0; i < fr->nre; i++)
466 dsqr(ee_sum[i].esum/nsum
467 - (ee_sum[i].esum + fr->ener[i].esum)/(nsum + fr->nsum))*
468 nsum*(nsum + fr->nsum)/(double)fr->nsum;
469 ee_sum[i].esum += fr->ener[i].esum;
472 nsteps += fr->nsteps;
479 fprintf(stderr, "\nWARNING: missing energy sums at time %f\n", fr->t);
485 *ee_sum_step = out_step;
486 *ee_sum_nsteps = nsteps;
490 int gmx_eneconv(int argc, char *argv[])
492 const char *desc[] = {
493 "With [IT]multiple files[it] specified for the [TT]-f[tt] option:[BR]",
494 "Concatenates several energy files in sorted order.",
495 "In the case of double time frames, the one",
496 "in the later file is used. By specifying [TT]-settime[tt] you will be",
497 "asked for the start time of each file. The input files are taken",
498 "from the command line,",
499 "such that the command [TT]eneconv -f *.edr -o fixed.edr[tt] should do",
501 "With [IT]one file[it] specified for [TT]-f[tt]:[BR]",
502 "Reads one energy file and writes another, applying the [TT]-dt[tt],",
503 "[TT]-offset[tt], [TT]-t0[tt] and [TT]-settime[tt] options and",
504 "converting to a different format if necessary (indicated by file",
506 "[TT]-settime[tt] is applied first, then [TT]-dt[tt]/[TT]-offset[tt]",
507 "followed by [TT]-b[tt] and [TT]-e[tt] to select which frames to write."
509 const char *bugs[] = {
510 "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."
512 ener_file_t in = NULL, out = NULL;
513 gmx_enxnm_t *enm = NULL;
515 ener_file_t in, out = NULL;
516 gmx_enxnm_t *enm = NULL;
518 t_enxframe *fr, *fro;
519 gmx_large_int_t ee_sum_step = 0, ee_sum_nsteps, ee_sum_nsum;
521 gmx_large_int_t lastfilestep, laststep, startstep, startstep_file = 0;
523 int nre, nremax, this_nre, nfile, f, i, j, kkk, nset, *set = NULL;
526 real *readtime, *settime, timestep, t1, tadjust;
527 char inputstring[STRLEN], *chptr, buf[22], buf2[22], buf3[22];
530 gmx_bool bNewFile, bFirst, bNewOutput;
532 gmx_bool warned_about_dh = FALSE;
533 t_enxblock *blocks = NULL;
535 int nblocks_alloc = 0;
538 { efEDR, "-f", NULL, ffRDMULT },
539 { efEDR, "-o", "fixed", ffWRITE },
542 #define NFILE asize(fnm)
544 static real delta_t = 0.0, toffset = 0, scalefac = 1;
545 static gmx_bool bSetTime = FALSE;
546 static gmx_bool bSort = TRUE, bError = TRUE;
547 static real begin = -1;
548 static real end = -1;
549 gmx_bool remove_dh = FALSE;
552 { "-b", FALSE, etREAL, {&begin},
553 "First time to use"},
554 { "-e", FALSE, etREAL, {&end},
556 { "-dt", FALSE, etREAL, {&delta_t},
557 "Only write out frame when t MOD dt = offset" },
558 { "-offset", FALSE, etREAL, {&toffset},
559 "Time offset for [TT]-dt[tt] option" },
560 { "-settime", FALSE, etBOOL, {&bSetTime},
561 "Change starting time interactively" },
562 { "-sort", FALSE, etBOOL, {&bSort},
563 "Sort energy files (not frames)"},
564 { "-rmdh", FALSE, etBOOL, {&remove_dh},
565 "Remove free energy block data" },
566 { "-scalefac", FALSE, etREAL, {&scalefac},
567 "Multiply energy component by this factor" },
568 { "-error", FALSE, etBOOL, {&bError},
569 "Stop on errors in the file" }
572 if (!parse_common_args(&argc, argv, PCA_BE_NICE, NFILE, fnm, asize(pa),
573 pa, asize(desc), desc, asize(bugs), bugs, &oenv))
584 laststep = startstep = 0;
586 nfile = opt2fns(&fnms, "-f", NFILE, fnm);
590 gmx_fatal(FARGS, "No input files!");
593 snew(settime, nfile+1);
594 snew(readtime, nfile+1);
595 snew(cont_type, nfile+1);
597 nre = scan_ene_files(fnms, nfile, readtime, ×tep, &nremax);
598 edit_files(fnms, nfile, readtime, settime, cont_type, bSetTime, bSort);
602 snew(ee_sum, nremax);
608 snew(fro->ener, nremax);
614 for (f = 0; f < nfile; f++)
618 in = open_enx(fnms[f], "r");
620 do_enxnms(in, &this_nre, &enm);
625 set = select_it(nre, enm, &nset);
628 /* write names to the output file */
629 out = open_enx(opt2fn("-o", NFILE, fnm), "w");
630 do_enxnms(out, &nre, &enm);
633 /* start reading from the next file */
634 while ((fro->t <= (settime[f+1] + GMX_REAL_EPS)) &&
639 startstep_file = fr->step;
640 tadjust = settime[f] - fr->t;
641 if (cont_type[f+1] == TIME_LAST)
643 settime[f+1] = readtime[f+1]-readtime[f]+settime[f];
644 cont_type[f+1] = TIME_EXPLICIT;
649 if (tadjust + fr->t <= last_t)
651 /* Skip this frame, since we already have it / past it */
654 fprintf(debug, "fr->step %s, fr->t %.4f\n",
655 gmx_step_str(fr->step, buf), fr->t);
656 fprintf(debug, "tadjust %12.6e + fr->t %12.6e <= t %12.6e\n",
657 tadjust, fr->t, last_t);
662 fro->step = lastfilestep + fr->step - startstep_file;
663 fro->t = tadjust + fr->t;
665 bWrite = ((begin < 0 || (begin >= 0 && (fro->t >= begin-GMX_REAL_EPS))) &&
666 (end < 0 || (end >= 0 && (fro->t <= end +GMX_REAL_EPS))) &&
667 (fro->t <= settime[f+1]+0.5*timestep));
672 "fr->step %s, fr->t %.4f, fro->step %s fro->t %.4f, w %d\n",
673 gmx_step_str(fr->step, buf), fr->t,
674 gmx_step_str(fro->step, buf2), fro->t, bWrite);
679 if ((end > 0) && (fro->t > end+GMX_REAL_EPS))
686 if (fro->t >= begin-GMX_REAL_EPS)
691 startstep = fr->step;
695 update_ee_sum(nre, &ee_sum_step, &ee_sum_nsteps, &ee_sum_nsum, ee_sum,
700 /* determine if we should write it */
701 if (bWrite && (delta_t == 0 || bRmod(fro->t, toffset, delta_t)))
703 laststep = fro->step;
708 fprintf(stderr, "\nContinue writing frames from t=%g, step=%s\n",
709 fro->t, gmx_step_str(fro->step, buf));
712 /* Copy the energies */
713 for (i = 0; i < nre; i++)
715 fro->ener[i].e = fr->ener[i].e;
718 fro->nsteps = ee_sum_nsteps;
721 if (ee_sum_nsum <= 1)
727 fro->nsum = gmx_large_int_to_int(ee_sum_nsum,
728 "energy average summation");
729 /* Copy the energy sums */
730 for (i = 0; i < nre; i++)
732 fro->ener[i].esum = ee_sum[i].esum;
733 fro->ener[i].eav = ee_sum[i].eav;
736 /* We wrote the energies, so reset the counts */
742 for (kkk = 0; kkk < nset; kkk++)
744 fro->ener[set[kkk]].e *= scalefac;
747 fro->ener[set[kkk]].eav *= scalefac*scalefac;
748 fro->ener[set[kkk]].esum *= scalefac;
752 /* Copy restraint stuff */
753 /*fro->ndisre = fr->ndisre;
754 fro->disre_rm3tav = fr->disre_rm3tav;
755 fro->disre_rt = fr->disre_rt;*/
756 fro->nblock = fr->nblock;
757 /*fro->nr = fr->nr;*/
758 fro->block = fr->block;
760 /* check if we have blocks with delta_h data and are throwing
767 if (!blocks || nblocks_alloc < fr->nblock)
769 /* we pre-allocate the blocks */
770 nblocks_alloc = fr->nblock;
771 snew(blocks, nblocks_alloc);
773 nblocks = 0; /* number of blocks so far */
775 for (i = 0; i < fr->nblock; i++)
777 if ( (fr->block[i].id != enxDHCOLL) &&
778 (fr->block[i].id != enxDH) &&
779 (fr->block[i].id != enxDHHIST) )
781 /* copy everything verbatim */
782 blocks[nblocks] = fr->block[i];
786 /* now set the block pointer to the new blocks */
787 fro->nblock = nblocks;
790 else if (delta_t > 0)
792 if (!warned_about_dh)
794 for (i = 0; i < fr->nblock; i++)
796 if (fr->block[i].id == enxDH ||
797 fr->block[i].id == enxDHHIST)
800 if (fr->block[i].id == enxDH)
802 size = fr->block[i].sub[2].nr;
810 printf("\nWARNING: %s contains delta H blocks or histograms for which\n"
811 " some data is thrown away on a block-by-block basis, where each block\n"
812 " contains up to %d samples.\n"
813 " This is almost certainly not what you want.\n"
814 " Use the -rmdh option to throw all delta H samples away.\n"
815 " Use g_energy -odh option to extract these samples.\n",
817 warned_about_dh = TRUE;
827 if (noutfr % 1000 == 0)
829 fprintf(stderr, "Writing frame time %g ", fro->t);
838 printf("\nLast step written from %s: t %g, step %s\n",
839 fnms[f], last_t, gmx_step_str(laststep, buf));
840 lastfilestep = laststep;
842 /* set the next time from the last in previous file */
843 if (cont_type[f+1] == TIME_CONTINUE)
845 settime[f+1] = fro->t;
846 /* in this case we have already written the last frame of
847 * previous file, so update begin to avoid doubling it
848 * with the start of the next file
850 begin = fro->t+0.5*timestep;
851 /* cont_type[f+1]==TIME_EXPLICIT; */
854 if ((fro->t < end) && (f < nfile-1) &&
855 (fro->t < settime[f+1]-1.5*timestep))
858 "\nWARNING: There might be a gap around t=%g\n", fro->t);
861 /* move energies to lastee */
863 free_enxnms(this_nre, enm);
865 fprintf(stderr, "\n");
869 fprintf(stderr, "No frames written.\n");
873 fprintf(stderr, "Last frame written was at step %s, time %f\n",
874 gmx_step_str(fro->step, buf), fro->t);
875 fprintf(stderr, "Wrote %d frames\n", noutfr);