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 "statutil.h"
  45 #include "disre.h"
  46 #include "names.h"
  47 #include "copyrite.h"
  48 #include "macros.h"
  49 #include "gmx_fatal.h"
  50 #include "enxio.h"
  51 #include "vec.h"
  52 #include "gmx_ana.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_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)
 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_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)-
 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_large_int_t *ee_sum_step,
 411                           gmx_large_int_t *ee_sum_nsteps,
 412                           gmx_large_int_t *ee_sum_nsum,
 413                           t_energy *ee_sum,
 414                           t_enxframe *fr, int out_step)
 415 {
 416     gmx_large_int_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]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_large_int_t ee_sum_step = 0, ee_sum_nsteps, ee_sum_nsum;
 521     t_energy       *ee_sum;
 522     gmx_large_int_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     parse_common_args(&argc, argv, PCA_BE_NICE, NFILE, fnm, asize(pa),
 574                       pa, asize(desc), desc, asize(bugs), bugs, &oenv);
 575     tadjust  = 0;
 576     nremax   = 0;
 577     nset     = 0;
 578     timestep = 0.0;
 579     snew(fnms, argc);
 580     nfile        = 0;
 581     lastfilestep = 0;
 582     laststep     = startstep = 0;
 583
 584     nfile = opt2fns(&fnms, "-f", NFILE, fnm);
 585
 586     if (!nfile)
 587     {
 588         gmx_fatal(FARGS, "No input files!");
 589     }
 590
 591     snew(settime, nfile+1);
 592     snew(readtime, nfile+1);
 593     snew(cont_type, nfile+1);
 594
 595     nre = scan_ene_files(fnms, nfile, readtime, &timestep, &nremax);
 596     edit_files(fnms, nfile, readtime, settime, cont_type, bSetTime, bSort);
 597
 598     ee_sum_nsteps = 0;
 599     ee_sum_nsum   = 0;
 600     snew(ee_sum, nremax);
 601
 602     snew(fr, 1);
 603     snew(fro, 1);
 604     fro->t   = -1e20;
 605     fro->nre = nre;
 606     snew(fro->ener, nremax);
 607
 608     noutfr = 0;
 609     bFirst = TRUE;
 610
 611     last_t = fro->t;
 612     for (f = 0; f < nfile; f++)
 613     {
 614         bNewFile   = TRUE;
 615         bNewOutput = TRUE;
 616         in         = open_enx(fnms[f], "r");
 617         enm        = NULL;
 618         do_enxnms(in, &this_nre, &enm);
 619         if (f == 0)
 620         {
 621             if (scalefac != 1)
 622             {
 623                 set = select_it(nre, enm, &nset);
 624             }
 625
 626             /* write names to the output file */
 627             out = open_enx(opt2fn("-o", NFILE, fnm), "w");
 628             do_enxnms(out, &nre, &enm);
 629         }
 630
 631         /* start reading from the next file */
 632         while ((fro->t <= (settime[f+1] + GMX_REAL_EPS)) &&
 633                do_enx(in, fr))
 634         {
 635             if (bNewFile)
 636             {
 637                 startstep_file = fr->step;
 638                 tadjust        = settime[f] - fr->t;
 639                 if (cont_type[f+1] == TIME_LAST)
 640                 {
 641                     settime[f+1]   = readtime[f+1]-readtime[f]+settime[f];
 642                     cont_type[f+1] = TIME_EXPLICIT;
 643                 }
 644                 bNewFile = FALSE;
 645             }
 646
 647             if (tadjust + fr->t <= last_t)
 648             {
 649                 /* Skip this frame, since we already have it / past it */
 650                 if (debug)
 651                 {
 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);
 656                 }
 657                 continue;
 658             }
 659
 660             fro->step = lastfilestep + fr->step - startstep_file;
 661             fro->t    = tadjust  + fr->t;
 662
 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));
 666
 667             if (debug)
 668             {
 669                 fprintf(debug,
 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);
 673             }
 674
 675             if (bError)
 676             {
 677                 if ((end > 0) && (fro->t > end+GMX_REAL_EPS))
 678                 {
 679                     f = nfile;
 680                     break;
 681                 }
 682             }
 683
 684             if (fro->t >= begin-GMX_REAL_EPS)
 685             {
 686                 if (bFirst)
 687                 {
 688                     bFirst    = FALSE;
 689                     startstep = fr->step;
 690                 }
 691                 if (bWrite)
 692                 {
 693                     update_ee_sum(nre, &ee_sum_step, &ee_sum_nsteps, &ee_sum_nsum, ee_sum,
 694                                   fr, fro->step);
 695                 }
 696             }
 697
 698             /* determine if we should write it */
 699             if (bWrite && (delta_t == 0 || bRmod(fro->t, toffset, delta_t)))
 700             {
 701                 laststep = fro->step;
 702                 last_t   = fro->t;
 703                 if (bNewOutput)
 704                 {
 705                     bNewOutput = FALSE;
 706                     fprintf(stderr, "\nContinue writing frames from t=%g, step=%s\n",
 707                             fro->t, gmx_step_str(fro->step, buf));
 708                 }
 709
 710                 /* Copy the energies */
 711                 for (i = 0; i < nre; i++)
 712                 {
 713                     fro->ener[i].e = fr->ener[i].e;
 714                 }
 715
 716                 fro->nsteps = ee_sum_nsteps;
 717                 fro->dt     = fr->dt;
 718
 719                 if (ee_sum_nsum <= 1)
 720                 {
 721                     fro->nsum = 0;
 722                 }
 723                 else
 724                 {
 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++)
 729                     {
 730                         fro->ener[i].esum = ee_sum[i].esum;
 731                         fro->ener[i].eav  = ee_sum[i].eav;
 732                     }
 733                 }
 734                 /* We wrote the energies, so reset the counts */
 735                 ee_sum_nsteps = 0;
 736                 ee_sum_nsum   = 0;
 737
 738                 if (scalefac != 1)
 739                 {
 740                     for (kkk = 0; kkk < nset; kkk++)
 741                     {
 742                         fro->ener[set[kkk]].e    *= scalefac;
 743                         if (fro->nsum > 0)
 744                         {
 745                             fro->ener[set[kkk]].eav  *= scalefac*scalefac;
 746                             fro->ener[set[kkk]].esum *= scalefac;
 747                         }
 748                     }
 749                 }
 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;
 757
 758                 /* check if we have blocks with delta_h data and are throwing
 759                    away data */
 760                 if (fro->nblock > 0)
 761                 {
 762                     if (remove_dh)
 763                     {
 764                         int i;
 765                         if (!blocks || nblocks_alloc < fr->nblock)
 766                         {
 767                             /* we pre-allocate the blocks */
 768                             nblocks_alloc = fr->nblock;
 769                             snew(blocks, nblocks_alloc);
 770                         }
 771                         nblocks = 0; /* number of blocks so far */
 772
 773                         for (i = 0; i < fr->nblock; i++)
 774                         {
 775                             if ( (fr->block[i].id != enxDHCOLL) &&
 776                                  (fr->block[i].id != enxDH) &&
 777                                  (fr->block[i].id != enxDHHIST) )
 778                             {
 779                                 /* copy everything verbatim */
 780                                 blocks[nblocks] = fr->block[i];
 781                                 nblocks++;
 782                             }
 783                         }
 784                         /* now set the block pointer to the new blocks */
 785                         fro->nblock = nblocks;
 786                         fro->block  = blocks;
 787                     }
 788                     else if (delta_t > 0)
 789                     {
 790                         if (!warned_about_dh)
 791                         {
 792                             for (i = 0; i < fr->nblock; i++)
 793                             {
 794                                 if (fr->block[i].id == enxDH ||
 795                                     fr->block[i].id == enxDHHIST)
 796                                 {
 797                                     int size;
 798                                     if (fr->block[i].id == enxDH)
 799                                     {
 800                                         size = fr->block[i].sub[2].nr;
 801                                     }
 802                                     else
 803                                     {
 804                                         size = fr->nsteps;
 805                                     }
 806                                     if (size > 0)
 807                                     {
 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",
 814                                                fnms[f], size);
 815                                         warned_about_dh = TRUE;
 816                                         break;
 817                                     }
 818                                 }
 819                             }
 820                         }
 821                     }
 822                 }
 823
 824                 do_enx(out, fro);
 825                 if (noutfr % 1000 == 0)
 826                 {
 827                     fprintf(stderr, "Writing frame time %g    ", fro->t);
 828                 }
 829                 noutfr++;
 830             }
 831         }
 832         if (f == nfile)
 833         {
 834             f--;
 835         }
 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;
 839
 840         /* set the next time from the last in previous file */
 841         if (cont_type[f+1] == TIME_CONTINUE)
 842         {
 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
 847              */
 848             begin = fro->t+0.5*timestep;
 849             /* cont_type[f+1]==TIME_EXPLICIT; */
 850         }
 851
 852         if ((fro->t < end) && (f < nfile-1) &&
 853             (fro->t < settime[f+1]-1.5*timestep))
 854         {
 855             fprintf(stderr,
 856                     "\nWARNING: There might be a gap around t=%g\n", fro->t);
 857         }
 858
 859         /* move energies to lastee */
 860         close_enx(in);
 861         free_enxnms(this_nre, enm);
 862
 863         fprintf(stderr, "\n");
 864     }
 865     if (noutfr == 0)
 866     {
 867         fprintf(stderr, "No frames written.\n");
 868     }
 869     else
 870     {
 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);
 874     }
 875
 876     thanx(stderr);
 877     return 0;
 878 }