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"
49 #include "gromacs/fileio/enxio.h"
52 #include "gromacs/fileio/trxio.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 if (!parse_common_args(&argc, argv, PCA_BE_NICE, NFILE, fnm, asize(pa),
574 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_large_int_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);