src/gromacs/gmxana/gmx_eneconv.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  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.
  14
  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.
  19  *
  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.
  26  *
  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.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Green Red Orange Magenta Azure Cyan Skyblue
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38 #include <string.h>
  39 #include <math.h>
  40
  41 #include "string2.h"
  42 #include "typedefs.h"
  43 #include "smalloc.h"
  44 #include "gromacs/commandline/pargs.h"
  45 #include "disre.h"
  46 #include "names.h"
  47 #include "macros.h"
  48 #include "gmx_fatal.h"
  49 #include "gromacs/fileio/enxio.h"
  50 #include "vec.h"
  51 #include "gmx_ana.h"
  52 #include "gromacs/fileio/trxio.h"
  53
  54 #define TIME_EXPLICIT 0
  55 #define TIME_CONTINUE 1
  56 #define TIME_LAST     2
  57 #ifndef FLT_MAX
  58 #define FLT_MAX 1e36
  59 #endif
  60
  61 static int *select_it(int nre, gmx_enxnm_t *nm, int *nset)
  62 {
  63     gmx_bool *bE;
  64     int       n, k, j, i;
  65     int      *set;
  66     gmx_bool  bVerbose = TRUE;
  67
  68     if ((getenv("VERBOSE")) != NULL)
  69     {
  70         bVerbose = FALSE;
  71     }
  72
  73     fprintf(stderr, "Select the terms you want to scale from the following list\n");
  74     fprintf(stderr, "End your selection with 0\n");
  75
  76     if (bVerbose)
  77     {
  78         for (k = 0; (k < nre); )
  79         {
  80             for (j = 0; (j < 4) && (k < nre); j++, k++)
  81             {
  82                 fprintf(stderr, " %3d=%14s", k+1, nm[k].name);
  83             }
  84             fprintf(stderr, "\n");
  85         }
  86     }
  87
  88     snew(bE, nre);
  89     do
  90     {
  91         if (1 != scanf("%d", &n))
  92         {
  93             gmx_fatal(FARGS, "Cannot read energy term");
  94         }
  95         if ((n > 0) && (n <= nre))
  96         {
  97             bE[n-1] = TRUE;
  98         }
  99     }
 100     while (n != 0);
 101
 102     snew(set, nre);
 103     for (i = (*nset) = 0; (i < nre); i++)
 104     {
 105         if (bE[i])
 106         {
 107             set[(*nset)++] = i;
 108         }
 109     }
 110
 111     sfree(bE);
 112
 113     return set;
 114 }
 115
 116 static void sort_files(char **fnms, real *settime, int nfile)
 117 {
 118     int   i, j, minidx;
 119     real  timeswap;
 120     char *chptr;
 121
 122     for (i = 0; i < nfile; i++)
 123     {
 124         minidx = i;
 125         for (j = i+1; j < nfile; j++)
 126         {
 127             if (settime[j] < settime[minidx])
 128             {
 129                 minidx = j;
 130             }
 131         }
 132         if (minidx != i)
 133         {
 134             timeswap        = settime[i];
 135             settime[i]      = settime[minidx];
 136             settime[minidx] = timeswap;
 137             chptr           = fnms[i];
 138             fnms[i]         = fnms[minidx];
 139             fnms[minidx]    = chptr;
 140         }
 141     }
 142 }
 143
 144
 145 static int scan_ene_files(char **fnms, int nfiles,
 146                           real *readtime, real *timestep, int *nremax)
 147 {
 148     /* Check number of energy terms and start time of all files */
 149     int          f, i, nre, nremin = 0, nresav = 0;
 150     ener_file_t  in;
 151     real         t1, t2;
 152     char         inputstring[STRLEN];
 153     gmx_enxnm_t *enm;
 154     t_enxframe  *fr;
 155
 156     snew(fr, 1);
 157
 158     for (f = 0; f < nfiles; f++)
 159     {
 160         in  = open_enx(fnms[f], "r");
 161         enm = NULL;
 162         do_enxnms(in, &nre, &enm);
 163
 164         if (f == 0)
 165         {
 166             nresav  = nre;
 167             nremin  = nre;
 168             *nremax = nre;
 169             do_enx(in, fr);
 170             t1 = fr->t;
 171             do_enx(in, fr);
 172             t2          = fr->t;
 173             *timestep   = t2-t1;
 174             readtime[f] = t1;
 175             close_enx(in);
 176         }
 177         else
 178         {
 179             nremin  = min(nremin, fr->nre);
 180             *nremax = max(*nremax, fr->nre);
 181             if (nre != nresav)
 182             {
 183                 fprintf(stderr,
 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);
 186                 fprintf(stderr,
 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))
 190                 {
 191                     gmx_fatal(FARGS, "Error reading user input");
 192                 }
 193                 if (inputstring[0] != 'y' && inputstring[0] != 'Y')
 194                 {
 195                     fprintf(stderr, "Will not convert\n");
 196                     exit(0);
 197                 }
 198                 nresav = fr->nre;
 199             }
 200             do_enx(in, fr);
 201             readtime[f] = fr->t;
 202             close_enx(in);
 203         }
 204         fprintf(stderr, "\n");
 205         free_enxnms(nre, enm);
 206     }
 207
 208     free_enxframe(fr);
 209     sfree(fr);
 210
 211     return nremin;
 212 }
 213
 214
 215 static void edit_files(char **fnms, int nfiles, real *readtime,
 216                        real *settime, int *cont_type, gmx_bool bSetTime, gmx_bool bSort)
 217 {
 218     int      i;
 219     gmx_bool ok;
 220     char     inputstring[STRLEN], *chptr;
 221
 222     if (bSetTime)
 223     {
 224         if (nfiles == 1)
 225         {
 226             fprintf(stderr, "\n\nEnter the new start time:\n\n");
 227         }
 228         else
 229         {
 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");
 239         }
 240
 241         fprintf(stderr, "          File             Current start       New start\n"
 242                 "---------------------------------------------------------\n");
 243
 244         for (i = 0; i < nfiles; i++)
 245         {
 246             fprintf(stderr, "%25s   %10.3f             ", fnms[i], readtime[i]);
 247             ok = FALSE;
 248             do
 249             {
 250                 if (NULL == fgets(inputstring, STRLEN-1, stdin))
 251                 {
 252                     gmx_fatal(FARGS, "Error reading user input");
 253                 }
 254                 inputstring[strlen(inputstring)-1] = 0;
 255
 256                 if (inputstring[0] == 'c' || inputstring[0] == 'C')
 257                 {
 258                     cont_type[i] = TIME_CONTINUE;
 259                     bSort        = FALSE;
 260                     ok           = TRUE;
 261                     settime[i]   = FLT_MAX;
 262                 }
 263                 else if (inputstring[0] == 'l' ||
 264                          inputstring[0] == 'L')
 265                 {
 266                     cont_type[i] = TIME_LAST;
 267                     bSort        = FALSE;
 268                     ok           = TRUE;
 269                     settime[i]   = FLT_MAX;
 270                 }
 271                 else
 272                 {
 273                     settime[i] = strtod(inputstring, &chptr);
 274                     if (chptr == inputstring)
 275                     {
 276                         fprintf(stderr, "Try that again: ");
 277                     }
 278                     else
 279                     {
 280                         cont_type[i] = TIME_EXPLICIT;
 281                         ok           = TRUE;
 282                     }
 283                 }
 284             }
 285             while (!ok);
 286         }
 287         if (cont_type[0] != TIME_EXPLICIT)
 288         {
 289             cont_type[0] = TIME_EXPLICIT;
 290             settime[0]   = 0;
 291         }
 292     }
 293     else
 294     {
 295         for (i = 0; i < nfiles; i++)
 296         {
 297             settime[i] = readtime[i];
 298         }
 299     }
 300
 301     if (bSort && (nfiles > 1))
 302     {
 303         sort_files(fnms, settime, nfiles);
 304     }
 305     else
 306     {
 307         fprintf(stderr, "Sorting disabled.\n");
 308     }
 309
 310
 311     /* Write out the new order and start times */
 312     fprintf(stderr, "\nSummary of files and start times used:\n\n"
 313             "          File                Start time\n"
 314             "-----------------------------------------\n");
 315     for (i = 0; i < nfiles; i++)
 316     {
 317         switch (cont_type[i])
 318         {
 319             case TIME_EXPLICIT:
 320                 fprintf(stderr, "%25s   %10.3f\n", fnms[i], settime[i]);
 321                 break;
 322             case TIME_CONTINUE:
 323                 fprintf(stderr, "%25s        Continue from end of last file\n", fnms[i]);
 324                 break;
 325             case TIME_LAST:
 326                 fprintf(stderr, "%25s        Change by same amount as last file\n", fnms[i]);
 327                 break;
 328         }
 329     }
 330     fprintf(stderr, "\n");
 331
 332     settime[nfiles]   = FLT_MAX;
 333     cont_type[nfiles] = TIME_EXPLICIT;
 334     readtime[nfiles]  = FLT_MAX;
 335 }
 336
 337
 338 static void copy_ee(t_energy *src, t_energy *dst, int nre)
 339 {
 340     int i;
 341
 342     for (i = 0; i < nre; i++)
 343     {
 344         dst[i].e    = src[i].e;
 345         dst[i].esum = src[i].esum;
 346         dst[i].eav  = src[i].eav;
 347     }
 348 }
 349
 350 static void update_ee(t_energy *lastee, gmx_int64_t laststep,
 351                       t_energy *startee, gmx_int64_t startstep,
 352                       t_energy *ee, int step,
 353                       t_energy *outee, int nre)
 354 {
 355     int    i;
 356     double sigmacorr, nom, denom;
 357     double prestart_esum;
 358     double prestart_sigma;
 359
 360     for (i = 0; i < nre; i++)
 361     {
 362         outee[i].e = ee[i].e;
 363         /* add statistics from earlier file if present */
 364         if (laststep > 0)
 365         {
 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;
 371         }
 372         else
 373         {
 374             /* otherwise just copy to output */
 375             outee[i].esum = ee[i].esum;
 376             outee[i].eav  = ee[i].eav;
 377         }
 378
 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.
 383          */
 384         if (startstep > 0)
 385         {
 386             gmx_int64_t q = laststep+step;
 387             gmx_int64_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)-
 393                 outee[i].esum/(q);
 394             outee[i].esum -= prestart_esum;
 395             if (q-p+1 > 0)
 396             {
 397                 outee[i].eav = outee[i].eav-prestart_sigma-
 398                     sigmacorr*sigmacorr*((p-1)*q)/(q-p+1);
 399             }
 400         }
 401
 402         if ((outee[i].eav/(laststep+step+1)) < (GMX_REAL_EPS))
 403         {
 404             outee[i].eav = 0;
 405         }
 406     }
 407 }
 408
 409 static void update_ee_sum(int nre,
 410                           gmx_int64_t *ee_sum_step,
 411                           gmx_int64_t *ee_sum_nsteps,
 412                           gmx_int64_t *ee_sum_nsum,
 413                           t_energy *ee_sum,
 414                           t_enxframe *fr, int out_step)
 415 {
 416     gmx_int64_t     nsteps, nsum, fr_nsum;
 417     int             i;
 418
 419     nsteps = *ee_sum_nsteps;
 420     nsum   = *ee_sum_nsum;
 421
 422     fr_nsum = fr->nsum;
 423     if (fr_nsum == 0)
 424     {
 425         fr_nsum = 1;
 426     }
 427
 428     if (nsteps == 0)
 429     {
 430         if (fr_nsum == 1)
 431         {
 432             for (i = 0; i < nre; i++)
 433             {
 434                 ee_sum[i].esum = fr->ener[i].e;
 435                 ee_sum[i].eav  = 0;
 436             }
 437         }
 438         else
 439         {
 440             for (i = 0; i < nre; i++)
 441             {
 442                 ee_sum[i].esum = fr->ener[i].esum;
 443                 ee_sum[i].eav  = fr->ener[i].eav;
 444             }
 445         }
 446         nsteps = fr->nsteps;
 447         nsum   = fr_nsum;
 448     }
 449     else if (out_step + *ee_sum_nsum - *ee_sum_step == nsteps + fr->nsteps)
 450     {
 451         if (fr_nsum == 1)
 452         {
 453             for (i = 0; i < nre; i++)
 454             {
 455                 ee_sum[i].eav  +=
 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;
 459             }
 460         }
 461         else
 462         {
 463             for (i = 0; i < fr->nre; i++)
 464             {
 465                 ee_sum[i].eav  +=
 466                     fr->ener[i].eav +
 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;
 471             }
 472         }
 473         nsteps += fr->nsteps;
 474         nsum   += fr_nsum;
 475     }
 476     else
 477     {
 478         if (fr->nsum != 0)
 479         {
 480             fprintf(stderr, "\nWARNING: missing energy sums at time %f\n", fr->t);
 481         }
 482         nsteps = 0;
 483         nsum   = 0;
 484     }
 485
 486     *ee_sum_step   = out_step;
 487     *ee_sum_nsteps = nsteps;
 488     *ee_sum_nsum   = nsum;
 489 }
 490
 491 int gmx_eneconv(int argc, char *argv[])
 492 {
 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]gmx eneconv -f *.edr -o fixed.edr[tt] should do",
 501         "the trick. [PAR]",
 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",
 506         "extentions).[PAR]",
 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."
 509     };
 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."
 512     };
 513     ener_file_t     in  = NULL, out = NULL;
 514     gmx_enxnm_t    *enm = NULL;
 515 #if 0
 516     ener_file_t     in, out = NULL;
 517     gmx_enxnm_t    *enm = NULL;
 518 #endif
 519     t_enxframe     *fr, *fro;
 520     gmx_int64_t     ee_sum_step = 0, ee_sum_nsteps, ee_sum_nsum;
 521     t_energy       *ee_sum;
 522     gmx_int64_t     lastfilestep, laststep, startstep, startstep_file = 0;
 523     int             noutfr;
 524     int             nre, nremax, this_nre, nfile, f, i, j, kkk, nset, *set = NULL;
 525     double          last_t;
 526     char          **fnms;
 527     real           *readtime, *settime, timestep, t1, tadjust;
 528     char            inputstring[STRLEN], *chptr, buf[22], buf2[22], buf3[22];
 529     gmx_bool        ok;
 530     int            *cont_type;
 531     gmx_bool        bNewFile, bFirst, bNewOutput;
 532     output_env_t    oenv;
 533     gmx_bool        warned_about_dh = FALSE;
 534     t_enxblock     *blocks          = NULL;
 535     int             nblocks         = 0;
 536     int             nblocks_alloc   = 0;
 537
 538     t_filenm        fnm[] = {
 539         { efEDR, "-f", NULL,    ffRDMULT },
 540         { efEDR, "-o", "fixed", ffWRITE  },
 541     };
 542
 543 #define NFILE asize(fnm)
 544     gmx_bool         bWrite;
 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;
 551
 552     t_pargs          pa[] = {
 553         { "-b",        FALSE, etREAL, {&begin},
 554           "First time to use"},
 555         { "-e",        FALSE, etREAL, {&end},
 556           "Last time to use"},
 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" }
 571     };
 572
 573     if (!parse_common_args(&argc, argv, PCA_BE_NICE, NFILE, fnm, asize(pa),
 574                            pa, asize(desc), desc, asize(bugs), bugs, &oenv))
 575     {
 576         return 0;
 577     }
 578     tadjust  = 0;
 579     nremax   = 0;
 580     nset     = 0;
 581     timestep = 0.0;
 582     snew(fnms, argc);
 583     nfile        = 0;
 584     lastfilestep = 0;
 585     laststep     = startstep = 0;
 586
 587     nfile = opt2fns(&fnms, "-f", NFILE, fnm);
 588
 589     if (!nfile)
 590     {
 591         gmx_fatal(FARGS, "No input files!");
 592     }
 593
 594     snew(settime, nfile+1);
 595     snew(readtime, nfile+1);
 596     snew(cont_type, nfile+1);
 597
 598     nre = scan_ene_files(fnms, nfile, readtime, &timestep, &nremax);
 599     edit_files(fnms, nfile, readtime, settime, cont_type, bSetTime, bSort);
 600
 601     ee_sum_nsteps = 0;
 602     ee_sum_nsum   = 0;
 603     snew(ee_sum, nremax);
 604
 605     snew(fr, 1);
 606     snew(fro, 1);
 607     fro->t   = -1e20;
 608     fro->nre = nre;
 609     snew(fro->ener, nremax);
 610
 611     noutfr = 0;
 612     bFirst = TRUE;
 613
 614     last_t = fro->t;
 615     for (f = 0; f < nfile; f++)
 616     {
 617         bNewFile   = TRUE;
 618         bNewOutput = TRUE;
 619         in         = open_enx(fnms[f], "r");
 620         enm        = NULL;
 621         do_enxnms(in, &this_nre, &enm);
 622         if (f == 0)
 623         {
 624             if (scalefac != 1)
 625             {
 626                 set = select_it(nre, enm, &nset);
 627             }
 628
 629             /* write names to the output file */
 630             out = open_enx(opt2fn("-o", NFILE, fnm), "w");
 631             do_enxnms(out, &nre, &enm);
 632         }
 633
 634         /* start reading from the next file */
 635         while ((fro->t <= (settime[f+1] + GMX_REAL_EPS)) &&
 636                do_enx(in, fr))
 637         {
 638             if (bNewFile)
 639             {
 640                 startstep_file = fr->step;
 641                 tadjust        = settime[f] - fr->t;
 642                 if (cont_type[f+1] == TIME_LAST)
 643                 {
 644                     settime[f+1]   = readtime[f+1]-readtime[f]+settime[f];
 645                     cont_type[f+1] = TIME_EXPLICIT;
 646                 }
 647                 bNewFile = FALSE;
 648             }
 649
 650             if (tadjust + fr->t <= last_t)
 651             {
 652                 /* Skip this frame, since we already have it / past it */
 653                 if (debug)
 654                 {
 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);
 659                 }
 660                 continue;
 661             }
 662
 663             fro->step = lastfilestep + fr->step - startstep_file;
 664             fro->t    = tadjust  + fr->t;
 665
 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));
 669
 670             if (debug)
 671             {
 672                 fprintf(debug,
 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);
 676             }
 677
 678             if (bError)
 679             {
 680                 if ((end > 0) && (fro->t > end+GMX_REAL_EPS))
 681                 {
 682                     f = nfile;
 683                     break;
 684                 }
 685             }
 686
 687             if (fro->t >= begin-GMX_REAL_EPS)
 688             {
 689                 if (bFirst)
 690                 {
 691                     bFirst    = FALSE;
 692                     startstep = fr->step;
 693                 }
 694                 if (bWrite)
 695                 {
 696                     update_ee_sum(nre, &ee_sum_step, &ee_sum_nsteps, &ee_sum_nsum, ee_sum,
 697                                   fr, fro->step);
 698                 }
 699             }
 700
 701             /* determine if we should write it */
 702             if (bWrite && (delta_t == 0 || bRmod(fro->t, toffset, delta_t)))
 703             {
 704                 laststep = fro->step;
 705                 last_t   = fro->t;
 706                 if (bNewOutput)
 707                 {
 708                     bNewOutput = FALSE;
 709                     fprintf(stderr, "\nContinue writing frames from t=%g, step=%s\n",
 710                             fro->t, gmx_step_str(fro->step, buf));
 711                 }
 712
 713                 /* Copy the energies */
 714                 for (i = 0; i < nre; i++)
 715                 {
 716                     fro->ener[i].e = fr->ener[i].e;
 717                 }
 718
 719                 fro->nsteps = ee_sum_nsteps;
 720                 fro->dt     = fr->dt;
 721
 722                 if (ee_sum_nsum <= 1)
 723                 {
 724                     fro->nsum = 0;
 725                 }
 726                 else
 727                 {
 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++)
 732                     {
 733                         fro->ener[i].esum = ee_sum[i].esum;
 734                         fro->ener[i].eav  = ee_sum[i].eav;
 735                     }
 736                 }
 737                 /* We wrote the energies, so reset the counts */
 738                 ee_sum_nsteps = 0;
 739                 ee_sum_nsum   = 0;
 740
 741                 if (scalefac != 1)
 742                 {
 743                     for (kkk = 0; kkk < nset; kkk++)
 744                     {
 745                         fro->ener[set[kkk]].e    *= scalefac;
 746                         if (fro->nsum > 0)
 747                         {
 748                             fro->ener[set[kkk]].eav  *= scalefac*scalefac;
 749                             fro->ener[set[kkk]].esum *= scalefac;
 750                         }
 751                     }
 752                 }
 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;
 760
 761                 /* check if we have blocks with delta_h data and are throwing
 762                    away data */
 763                 if (fro->nblock > 0)
 764                 {
 765                     if (remove_dh)
 766                     {
 767                         int i;
 768                         if (!blocks || nblocks_alloc < fr->nblock)
 769                         {
 770                             /* we pre-allocate the blocks */
 771                             nblocks_alloc = fr->nblock;
 772                             snew(blocks, nblocks_alloc);
 773                         }
 774                         nblocks = 0; /* number of blocks so far */
 775
 776                         for (i = 0; i < fr->nblock; i++)
 777                         {
 778                             if ( (fr->block[i].id != enxDHCOLL) &&
 779                                  (fr->block[i].id != enxDH) &&
 780                                  (fr->block[i].id != enxDHHIST) )
 781                             {
 782                                 /* copy everything verbatim */
 783                                 blocks[nblocks] = fr->block[i];
 784                                 nblocks++;
 785                             }
 786                         }
 787                         /* now set the block pointer to the new blocks */
 788                         fro->nblock = nblocks;
 789                         fro->block  = blocks;
 790                     }
 791                     else if (delta_t > 0)
 792                     {
 793                         if (!warned_about_dh)
 794                         {
 795                             for (i = 0; i < fr->nblock; i++)
 796                             {
 797                                 if (fr->block[i].id == enxDH ||
 798                                     fr->block[i].id == enxDHHIST)
 799                                 {
 800                                     int size;
 801                                     if (fr->block[i].id == enxDH)
 802                                     {
 803                                         size = fr->block[i].sub[2].nr;
 804                                     }
 805                                     else
 806                                     {
 807                                         size = fr->nsteps;
 808                                     }
 809                                     if (size > 0)
 810                                     {
 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",
 817                                                fnms[f], size);
 818                                         warned_about_dh = TRUE;
 819                                         break;
 820                                     }
 821                                 }
 822                             }
 823                         }
 824                     }
 825                 }
 826
 827                 do_enx(out, fro);
 828                 if (noutfr % 1000 == 0)
 829                 {
 830                     fprintf(stderr, "Writing frame time %g    ", fro->t);
 831                 }
 832                 noutfr++;
 833             }
 834         }
 835         if (f == nfile)
 836         {
 837             f--;
 838         }
 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;
 842
 843         /* set the next time from the last in previous file */
 844         if (cont_type[f+1] == TIME_CONTINUE)
 845         {
 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
 850              */
 851             begin = fro->t+0.5*timestep;
 852             /* cont_type[f+1]==TIME_EXPLICIT; */
 853         }
 854
 855         if ((fro->t < end) && (f < nfile-1) &&
 856             (fro->t < settime[f+1]-1.5*timestep))
 857         {
 858             fprintf(stderr,
 859                     "\nWARNING: There might be a gap around t=%g\n", fro->t);
 860         }
 861
 862         /* move energies to lastee */
 863         close_enx(in);
 864         free_enxnms(this_nre, enm);
 865
 866         fprintf(stderr, "\n");
 867     }
 868     if (noutfr == 0)
 869     {
 870         fprintf(stderr, "No frames written.\n");
 871     }
 872     else
 873     {
 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);
 877     }
 878
 879     return 0;
 880 }