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
49 #include "gmx_fatal.h"
54 #define TIME_EXPLICIT 0
55 #define TIME_CONTINUE 1
61 static int *select_it(int nre, gmx_enxnm_t *nm, int *nset)
66 gmx_bool bVerbose = TRUE;
68 if ((getenv("VERBOSE")) != NULL)
73 fprintf(stderr, "Select the terms you want to scale from the following list\n");
74 fprintf(stderr, "End your selection with 0\n");
78 for (k = 0; (k < nre); )
80 for (j = 0; (j < 4) && (k < nre); j++, k++)
82 fprintf(stderr, " %3d=%14s", k+1, nm[k].name);
84 fprintf(stderr, "\n");
91 if (1 != scanf("%d", &n))
93 gmx_fatal(FARGS, "Cannot read energy term");
95 if ((n > 0) && (n <= nre))
103 for (i = (*nset) = 0; (i < nre); i++)
116 static void sort_files(char **fnms, real *settime, int nfile)
122 for (i = 0; i < nfile; i++)
125 for (j = i+1; j < nfile; j++)
127 if (settime[j] < settime[minidx])
134 timeswap = settime[i];
135 settime[i] = settime[minidx];
136 settime[minidx] = timeswap;
138 fnms[i] = fnms[minidx];
139 fnms[minidx] = chptr;
145 static int scan_ene_files(char **fnms, int nfiles,
146 real *readtime, real *timestep, int *nremax)
148 /* Check number of energy terms and start time of all files */
149 int f, i, nre, nremin = 0, nresav = 0;
152 char inputstring[STRLEN];
158 for (f = 0; f < nfiles; f++)
160 in = open_enx(fnms[f], "r");
162 do_enxnms(in, &nre, &enm);
179 nremin = min(nremin, fr->nre);
180 *nremax = max(*nremax, fr->nre);
184 "Energy files don't match, different number of energies:\n"
185 " %s: %d\n %s: %d\n", fnms[f-1], nresav, fnms[f], fr->nre);
187 "\nContinue conversion using only the first %d terms (n/y)?\n"
188 "(you should be sure that the energy terms match)\n", nremin);
189 if (NULL == fgets(inputstring, STRLEN-1, stdin))
191 gmx_fatal(FARGS, "Error reading user input");
193 if (inputstring[0] != 'y' && inputstring[0] != 'Y')
195 fprintf(stderr, "Will not convert\n");
204 fprintf(stderr, "\n");
205 free_enxnms(nre, enm);
215 static void edit_files(char **fnms, int nfiles, real *readtime,
216 real *settime, int *cont_type, gmx_bool bSetTime, gmx_bool bSort)
220 char inputstring[STRLEN], *chptr;
226 fprintf(stderr, "\n\nEnter the new start time:\n\n");
230 fprintf(stderr, "\n\nEnter the new start time for each file.\n"
231 "There are two special options, both disables sorting:\n\n"
232 "c (continue) - The start time is taken from the end\n"
233 "of the previous file. Use it when your continuation run\n"
234 "restarts with t=0 and there is no overlap.\n\n"
235 "l (last) - The time in this file will be changed the\n"
236 "same amount as in the previous. Use it when the time in the\n"
237 "new run continues from the end of the previous one,\n"
238 "since this takes possible overlap into account.\n\n");
241 fprintf(stderr, " File Current start New start\n"
242 "---------------------------------------------------------\n");
244 for (i = 0; i < nfiles; i++)
246 fprintf(stderr, "%25s %10.3f ", fnms[i], readtime[i]);
250 if (NULL == fgets(inputstring, STRLEN-1, stdin))
252 gmx_fatal(FARGS, "Error reading user input");
254 inputstring[strlen(inputstring)-1] = 0;
256 if (inputstring[0] == 'c' || inputstring[0] == 'C')
258 cont_type[i] = TIME_CONTINUE;
261 settime[i] = FLT_MAX;
263 else if (inputstring[0] == 'l' ||
264 inputstring[0] == 'L')
266 cont_type[i] = TIME_LAST;
269 settime[i] = FLT_MAX;
273 settime[i] = strtod(inputstring, &chptr);
274 if (chptr == inputstring)
276 fprintf(stderr, "Try that again: ");
280 cont_type[i] = TIME_EXPLICIT;
287 if (cont_type[0] != TIME_EXPLICIT)
289 cont_type[0] = TIME_EXPLICIT;
295 for (i = 0; i < nfiles; i++)
297 settime[i] = readtime[i];
301 if (bSort && (nfiles > 1))
303 sort_files(fnms, settime, nfiles);
307 fprintf(stderr, "Sorting disabled.\n");
311 /* Write out the new order and start times */
312 fprintf(stderr, "\nSummary of files and start times used:\n\n"
314 "-----------------------------------------\n");
315 for (i = 0; i < nfiles; i++)
317 switch (cont_type[i])
320 fprintf(stderr, "%25s %10.3f\n", fnms[i], settime[i]);
323 fprintf(stderr, "%25s Continue from end of last file\n", fnms[i]);
326 fprintf(stderr, "%25s Change by same amount as last file\n", fnms[i]);
330 fprintf(stderr, "\n");
332 settime[nfiles] = FLT_MAX;
333 cont_type[nfiles] = TIME_EXPLICIT;
334 readtime[nfiles] = FLT_MAX;
338 static void copy_ee(t_energy *src, t_energy *dst, int nre)
342 for (i = 0; i < nre; i++)
345 dst[i].esum = src[i].esum;
346 dst[i].eav = src[i].eav;
350 static void update_ee(t_energy *lastee, gmx_large_int_t laststep,
351 t_energy *startee, gmx_large_int_t startstep,
352 t_energy *ee, int step,
353 t_energy *outee, int nre)
356 double sigmacorr, nom, denom;
357 double prestart_esum;
358 double prestart_sigma;
360 for (i = 0; i < nre; i++)
362 outee[i].e = ee[i].e;
363 /* add statistics from earlier file if present */
366 outee[i].esum = lastee[i].esum+ee[i].esum;
367 nom = (lastee[i].esum*(step+1)-ee[i].esum*(laststep));
368 denom = ((step+1.0)*(laststep)*(step+1.0+laststep));
369 sigmacorr = nom*nom/denom;
370 outee[i].eav = lastee[i].eav+ee[i].eav+sigmacorr;
374 /* otherwise just copy to output */
375 outee[i].esum = ee[i].esum;
376 outee[i].eav = ee[i].eav;
379 /* if we didnt start to write at the first frame
380 * we must compensate the statistics for this
381 * there are laststep frames in the earlier file
382 * and step+1 frames in this one.
386 gmx_large_int_t q = laststep+step;
387 gmx_large_int_t p = startstep+1;
388 prestart_esum = startee[i].esum-startee[i].e;
389 sigmacorr = prestart_esum-(p-1)*startee[i].e;
390 prestart_sigma = startee[i].eav-
391 sigmacorr*sigmacorr/(p*(p-1));
392 sigmacorr = prestart_esum/(p-1)-
394 outee[i].esum -= prestart_esum;
397 outee[i].eav = outee[i].eav-prestart_sigma-
398 sigmacorr*sigmacorr*((p-1)*q)/(q-p+1);
402 if ((outee[i].eav/(laststep+step+1)) < (GMX_REAL_EPS))
409 static void update_ee_sum(int nre,
410 gmx_large_int_t *ee_sum_step,
411 gmx_large_int_t *ee_sum_nsteps,
412 gmx_large_int_t *ee_sum_nsum,
414 t_enxframe *fr, int out_step)
416 gmx_large_int_t nsteps, nsum, fr_nsum;
419 nsteps = *ee_sum_nsteps;
432 for (i = 0; i < nre; i++)
434 ee_sum[i].esum = fr->ener[i].e;
440 for (i = 0; i < nre; i++)
442 ee_sum[i].esum = fr->ener[i].esum;
443 ee_sum[i].eav = fr->ener[i].eav;
449 else if (out_step + *ee_sum_nsum - *ee_sum_step == nsteps + fr->nsteps)
453 for (i = 0; i < nre; i++)
456 dsqr(ee_sum[i].esum/nsum
457 - (ee_sum[i].esum + fr->ener[i].e)/(nsum + 1))*nsum*(nsum + 1);
458 ee_sum[i].esum += fr->ener[i].e;
463 for (i = 0; i < fr->nre; i++)
467 dsqr(ee_sum[i].esum/nsum
468 - (ee_sum[i].esum + fr->ener[i].esum)/(nsum + fr->nsum))*
469 nsum*(nsum + fr->nsum)/(double)fr->nsum;
470 ee_sum[i].esum += fr->ener[i].esum;
473 nsteps += fr->nsteps;
480 fprintf(stderr, "\nWARNING: missing energy sums at time %f\n", fr->t);
486 *ee_sum_step = out_step;
487 *ee_sum_nsteps = nsteps;
491 int gmx_eneconv(int argc, char *argv[])
493 const char *desc[] = {
494 "With [IT]multiple files[it] specified for the [TT]-f[tt] option:[BR]",
495 "Concatenates several energy files in sorted order.",
496 "In the case of double time frames, the one",
497 "in the later file is used. By specifying [TT]-settime[tt] you will be",
498 "asked for the start time of each file. The input files are taken",
499 "from the command line,",
500 "such that the command [TT]eneconv -f *.edr -o fixed.edr[tt] should do",
502 "With [IT]one file[it] specified for [TT]-f[tt]:[BR]",
503 "Reads one energy file and writes another, applying the [TT]-dt[tt],",
504 "[TT]-offset[tt], [TT]-t0[tt] and [TT]-settime[tt] options and",
505 "converting to a different format if necessary (indicated by file",
507 "[TT]-settime[tt] is applied first, then [TT]-dt[tt]/[TT]-offset[tt]",
508 "followed by [TT]-b[tt] and [TT]-e[tt] to select which frames to write."
510 const char *bugs[] = {
511 "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."
513 ener_file_t in = NULL, out = NULL;
514 gmx_enxnm_t *enm = NULL;
516 ener_file_t in, out = NULL;
517 gmx_enxnm_t *enm = NULL;
519 t_enxframe *fr, *fro;
520 gmx_large_int_t ee_sum_step = 0, ee_sum_nsteps, ee_sum_nsum;
522 gmx_large_int_t lastfilestep, laststep, startstep, startstep_file = 0;
524 int nre, nremax, this_nre, nfile, f, i, j, kkk, nset, *set = NULL;
527 real *readtime, *settime, timestep, t1, tadjust;
528 char inputstring[STRLEN], *chptr, buf[22], buf2[22], buf3[22];
531 gmx_bool bNewFile, bFirst, bNewOutput;
533 gmx_bool warned_about_dh = FALSE;
534 t_enxblock *blocks = NULL;
536 int nblocks_alloc = 0;
539 { efEDR, "-f", NULL, ffRDMULT },
540 { efEDR, "-o", "fixed", ffWRITE },
543 #define NFILE asize(fnm)
545 static real delta_t = 0.0, toffset = 0, scalefac = 1;
546 static gmx_bool bSetTime = FALSE;
547 static gmx_bool bSort = TRUE, bError = TRUE;
548 static real begin = -1;
549 static real end = -1;
550 gmx_bool remove_dh = FALSE;
553 { "-b", FALSE, etREAL, {&begin},
554 "First time to use"},
555 { "-e", FALSE, etREAL, {&end},
557 { "-dt", FALSE, etREAL, {&delta_t},
558 "Only write out frame when t MOD dt = offset" },
559 { "-offset", FALSE, etREAL, {&toffset},
560 "Time offset for [TT]-dt[tt] option" },
561 { "-settime", FALSE, etBOOL, {&bSetTime},
562 "Change starting time interactively" },
563 { "-sort", FALSE, etBOOL, {&bSort},
564 "Sort energy files (not frames)"},
565 { "-rmdh", FALSE, etBOOL, {&remove_dh},
566 "Remove free energy block data" },
567 { "-scalefac", FALSE, etREAL, {&scalefac},
568 "Multiply energy component by this factor" },
569 { "-error", FALSE, etBOOL, {&bError},
570 "Stop on errors in the file" }
573 parse_common_args(&argc, argv, PCA_BE_NICE, NFILE, fnm, asize(pa),
574 pa, asize(desc), desc, asize(bugs), bugs, &oenv);
582 laststep = startstep = 0;
584 nfile = opt2fns(&fnms, "-f", NFILE, fnm);
588 gmx_fatal(FARGS, "No input files!");
591 snew(settime, nfile+1);
592 snew(readtime, nfile+1);
593 snew(cont_type, nfile+1);
595 nre = scan_ene_files(fnms, nfile, readtime, ×tep, &nremax);
596 edit_files(fnms, nfile, readtime, settime, cont_type, bSetTime, bSort);
600 snew(ee_sum, nremax);
606 snew(fro->ener, nremax);
612 for (f = 0; f < nfile; f++)
616 in = open_enx(fnms[f], "r");
618 do_enxnms(in, &this_nre, &enm);
623 set = select_it(nre, enm, &nset);
626 /* write names to the output file */
627 out = open_enx(opt2fn("-o", NFILE, fnm), "w");
628 do_enxnms(out, &nre, &enm);
631 /* start reading from the next file */
632 while ((fro->t <= (settime[f+1] + GMX_REAL_EPS)) &&
637 startstep_file = fr->step;
638 tadjust = settime[f] - fr->t;
639 if (cont_type[f+1] == TIME_LAST)
641 settime[f+1] = readtime[f+1]-readtime[f]+settime[f];
642 cont_type[f+1] = TIME_EXPLICIT;
647 if (tadjust + fr->t <= last_t)
649 /* Skip this frame, since we already have it / past it */
652 fprintf(debug, "fr->step %s, fr->t %.4f\n",
653 gmx_step_str(fr->step, buf), fr->t);
654 fprintf(debug, "tadjust %12.6e + fr->t %12.6e <= t %12.6e\n",
655 tadjust, fr->t, last_t);
660 fro->step = lastfilestep + fr->step - startstep_file;
661 fro->t = tadjust + fr->t;
663 bWrite = ((begin < 0 || (begin >= 0 && (fro->t >= begin-GMX_REAL_EPS))) &&
664 (end < 0 || (end >= 0 && (fro->t <= end +GMX_REAL_EPS))) &&
665 (fro->t <= settime[f+1]+0.5*timestep));
670 "fr->step %s, fr->t %.4f, fro->step %s fro->t %.4f, w %d\n",
671 gmx_step_str(fr->step, buf), fr->t,
672 gmx_step_str(fro->step, buf2), fro->t, bWrite);
677 if ((end > 0) && (fro->t > end+GMX_REAL_EPS))
684 if (fro->t >= begin-GMX_REAL_EPS)
689 startstep = fr->step;
693 update_ee_sum(nre, &ee_sum_step, &ee_sum_nsteps, &ee_sum_nsum, ee_sum,
698 /* determine if we should write it */
699 if (bWrite && (delta_t == 0 || bRmod(fro->t, toffset, delta_t)))
701 laststep = fro->step;
706 fprintf(stderr, "\nContinue writing frames from t=%g, step=%s\n",
707 fro->t, gmx_step_str(fro->step, buf));
710 /* Copy the energies */
711 for (i = 0; i < nre; i++)
713 fro->ener[i].e = fr->ener[i].e;
716 fro->nsteps = ee_sum_nsteps;
719 if (ee_sum_nsum <= 1)
725 fro->nsum = gmx_large_int_to_int(ee_sum_nsum,
726 "energy average summation");
727 /* Copy the energy sums */
728 for (i = 0; i < nre; i++)
730 fro->ener[i].esum = ee_sum[i].esum;
731 fro->ener[i].eav = ee_sum[i].eav;
734 /* We wrote the energies, so reset the counts */
740 for (kkk = 0; kkk < nset; kkk++)
742 fro->ener[set[kkk]].e *= scalefac;
745 fro->ener[set[kkk]].eav *= scalefac*scalefac;
746 fro->ener[set[kkk]].esum *= scalefac;
750 /* Copy restraint stuff */
751 /*fro->ndisre = fr->ndisre;
752 fro->disre_rm3tav = fr->disre_rm3tav;
753 fro->disre_rt = fr->disre_rt;*/
754 fro->nblock = fr->nblock;
755 /*fro->nr = fr->nr;*/
756 fro->block = fr->block;
758 /* check if we have blocks with delta_h data and are throwing
765 if (!blocks || nblocks_alloc < fr->nblock)
767 /* we pre-allocate the blocks */
768 nblocks_alloc = fr->nblock;
769 snew(blocks, nblocks_alloc);
771 nblocks = 0; /* number of blocks so far */
773 for (i = 0; i < fr->nblock; i++)
775 if ( (fr->block[i].id != enxDHCOLL) &&
776 (fr->block[i].id != enxDH) &&
777 (fr->block[i].id != enxDHHIST) )
779 /* copy everything verbatim */
780 blocks[nblocks] = fr->block[i];
784 /* now set the block pointer to the new blocks */
785 fro->nblock = nblocks;
788 else if (delta_t > 0)
790 if (!warned_about_dh)
792 for (i = 0; i < fr->nblock; i++)
794 if (fr->block[i].id == enxDH ||
795 fr->block[i].id == enxDHHIST)
798 if (fr->block[i].id == enxDH)
800 size = fr->block[i].sub[2].nr;
808 printf("\nWARNING: %s contains delta H blocks or histograms for which\n"
809 " some data is thrown away on a block-by-block basis, where each block\n"
810 " contains up to %d samples.\n"
811 " This is almost certainly not what you want.\n"
812 " Use the -rmdh option to throw all delta H samples away.\n"
813 " Use g_energy -odh option to extract these samples.\n",
815 warned_about_dh = TRUE;
825 if (noutfr % 1000 == 0)
827 fprintf(stderr, "Writing frame time %g ", fro->t);
836 printf("\nLast step written from %s: t %g, step %s\n",
837 fnms[f], last_t, gmx_step_str(laststep, buf));
838 lastfilestep = laststep;
840 /* set the next time from the last in previous file */
841 if (cont_type[f+1] == TIME_CONTINUE)
843 settime[f+1] = fro->t;
844 /* in this case we have already written the last frame of
845 * previous file, so update begin to avoid doubling it
846 * with the start of the next file
848 begin = fro->t+0.5*timestep;
849 /* cont_type[f+1]==TIME_EXPLICIT; */
852 if ((fro->t < end) && (f < nfile-1) &&
853 (fro->t < settime[f+1]-1.5*timestep))
856 "\nWARNING: There might be a gap around t=%g\n", fro->t);
859 /* move energies to lastee */
861 free_enxnms(this_nre, enm);
863 fprintf(stderr, "\n");
867 fprintf(stderr, "No frames written.\n");
871 fprintf(stderr, "Last frame written was at step %s, time %f\n",
872 gmx_step_str(fro->step, buf), fro->t);
873 fprintf(stderr, "Wrote %d frames\n", noutfr);