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 "macros.h"
  48 #include "gmx_fatal.h"
  49 #include "enxio.h"
  50 #include "vec.h"
  51 #include "gmx_ana.h"
  52
  53 #define TIME_EXPLICIT 0
  54 #define TIME_CONTINUE 1
  55 #define TIME_LAST     2
  56 #ifndef FLT_MAX
  57 #define FLT_MAX 1e36
  58 #endif
  59
  60 static int *select_it(int nre, gmx_enxnm_t *nm, int *nset)
  61 {
  62     gmx_bool *bE;
  63     int       n, k, j, i;
  64     int      *set;
  65     gmx_bool  bVerbose = TRUE;
  66
  67     if ((getenv("VERBOSE")) != NULL)
  68     {
  69         bVerbose = FALSE;
  70     }
  71
  72     fprintf(stderr, "Select the terms you want to scale from the following list\n");
  73     fprintf(stderr, "End your selection with 0\n");
  74
  75     if (bVerbose)
  76     {
  77         for (k = 0; (k < nre); )
  78         {
  79             for (j = 0; (j < 4) && (k < nre); j++, k++)
  80             {
  81                 fprintf(stderr, " %3d=%14s", k+1, nm[k].name);
  82             }
  83             fprintf(stderr, "\n");
  84         }
  85     }
  86
  87     snew(bE, nre);
  88     do
  89     {
  90         if (1 != scanf("%d", &n))
  91         {
  92             gmx_fatal(FARGS, "Cannot read energy term");
  93         }
  94         if ((n > 0) && (n <= nre))
  95         {
  96             bE[n-1] = TRUE;
  97         }
  98     }
  99     while (n != 0);
 100
 101     snew(set, nre);
 102     for (i = (*nset) = 0; (i < nre); i++)
 103     {
 104         if (bE[i])
 105         {
 106             set[(*nset)++] = i;
 107         }
 108     }
 109
 110     sfree(bE);
 111
 112     return set;
 113 }
 114
 115 static void sort_files(char **fnms, real *settime, int nfile)
 116 {
 117     int   i, j, minidx;
 118     real  timeswap;
 119     char *chptr;
 120
 121     for (i = 0; i < nfile; i++)
 122     {
 123         minidx = i;
 124         for (j = i+1; j < nfile; j++)
 125         {
 126             if (settime[j] < settime[minidx])
 127             {
 128                 minidx = j;
 129             }
 130         }
 131         if (minidx != i)
 132         {
 133             timeswap        = settime[i];
 134             settime[i]      = settime[minidx];
 135             settime[minidx] = timeswap;
 136             chptr           = fnms[i];
 137             fnms[i]         = fnms[minidx];
 138             fnms[minidx]    = chptr;
 139         }
 140     }
 141 }
 142
 143
 144 static int scan_ene_files(char **fnms, int nfiles,
 145                           real *readtime, real *timestep, int *nremax)
 146 {
 147     /* Check number of energy terms and start time of all files */
 148     int          f, i, nre, nremin = 0, nresav = 0;
 149     ener_file_t  in;
 150     real         t1, t2;
 151     char         inputstring[STRLEN];
 152     gmx_enxnm_t *enm;
 153     t_enxframe  *fr;
 154
 155     snew(fr, 1);
 156
 157     for (f = 0; f < nfiles; f++)
 158     {
 159         in  = open_enx(fnms[f], "r");
 160         enm = NULL;
 161         do_enxnms(in, &nre, &enm);
 162
 163         if (f == 0)
 164         {
 165             nresav  = nre;
 166             nremin  = nre;
 167             *nremax = nre;
 168             do_enx(in, fr);
 169             t1 = fr->t;
 170             do_enx(in, fr);
 171             t2          = fr->t;
 172             *timestep   = t2-t1;
 173             readtime[f] = t1;
 174             close_enx(in);
 175         }
 176         else
 177         {
 178             nremin  = min(nremin, fr->nre);
 179             *nremax = max(*nremax, fr->nre);
 180             if (nre != nresav)
 181             {
 182                 fprintf(stderr,
 183                         "Energy files don't match, different number of energies:\n"
 184                         " %s: %d\n %s: %d\n", fnms[f-1], nresav, fnms[f], fr->nre);
 185                 fprintf(stderr,
 186                         "\nContinue conversion using only the first %d terms (n/y)?\n"
 187                         "(you should be sure that the energy terms match)\n", nremin);
 188                 if (NULL == fgets(inputstring, STRLEN-1, stdin))
 189                 {
 190                     gmx_fatal(FARGS, "Error reading user input");
 191                 }
 192                 if (inputstring[0] != 'y' && inputstring[0] != 'Y')
 193                 {
 194                     fprintf(stderr, "Will not convert\n");
 195                     exit(0);
 196                 }
 197                 nresav = fr->nre;
 198             }
 199             do_enx(in, fr);
 200             readtime[f] = fr->t;
 201             close_enx(in);
 202         }
 203         fprintf(stderr, "\n");
 204         free_enxnms(nre, enm);
 205     }
 206
 207     free_enxframe(fr);
 208     sfree(fr);
 209
 210     return nremin;
 211 }
 212
 213
 214 static void edit_files(char **fnms, int nfiles, real *readtime,
 215                        real *settime, int *cont_type, gmx_bool bSetTime, gmx_bool bSort)
 216 {
 217     int      i;
 218     gmx_bool ok;
 219     char     inputstring[STRLEN], *chptr;
 220
 221     if (bSetTime)
 222     {
 223         if (nfiles == 1)
 224         {
 225             fprintf(stderr, "\n\nEnter the new start time:\n\n");
 226         }
 227         else
 228         {
 229             fprintf(stderr, "\n\nEnter the new start time for each file.\n"
 230                     "There are two special options, both disables sorting:\n\n"
 231                     "c (continue) - The start time is taken from the end\n"
 232                     "of the previous file. Use it when your continuation run\n"
 233                     "restarts with t=0 and there is no overlap.\n\n"
 234                     "l (last) - The time in this file will be changed the\n"
 235                     "same amount as in the previous. Use it when the time in the\n"
 236                     "new run continues from the end of the previous one,\n"
 237                     "since this takes possible overlap into account.\n\n");
 238         }
 239
 240         fprintf(stderr, "          File             Current start       New start\n"
 241                 "---------------------------------------------------------\n");
 242
 243         for (i = 0; i < nfiles; i++)
 244         {
 245             fprintf(stderr, "%25s   %10.3f             ", fnms[i], readtime[i]);
 246             ok = FALSE;
 247             do
 248             {
 249                 if (NULL == fgets(inputstring, STRLEN-1, stdin))
 250                 {
 251                     gmx_fatal(FARGS, "Error reading user input");
 252                 }
 253                 inputstring[strlen(inputstring)-1] = 0;
 254
 255                 if (inputstring[0] == 'c' || inputstring[0] == 'C')
 256                 {
 257                     cont_type[i] = TIME_CONTINUE;
 258                     bSort        = FALSE;
 259                     ok           = TRUE;
 260                     settime[i]   = FLT_MAX;
 261                 }
 262                 else if (inputstring[0] == 'l' ||
 263                          inputstring[0] == 'L')
 264                 {
 265                     cont_type[i] = TIME_LAST;
 266                     bSort        = FALSE;
 267                     ok           = TRUE;
 268                     settime[i]   = FLT_MAX;
 269                 }
 270                 else
 271                 {
 272                     settime[i] = strtod(inputstring, &chptr);
 273                     if (chptr == inputstring)
 274                     {
 275                         fprintf(stderr, "Try that again: ");
 276                     }
 277                     else
 278                     {
 279                         cont_type[i] = TIME_EXPLICIT;
 280                         ok           = TRUE;
 281                     }
 282                 }
 283             }
 284             while (!ok);
 285         }
 286         if (cont_type[0] != TIME_EXPLICIT)
 287         {
 288             cont_type[0] = TIME_EXPLICIT;
 289             settime[0]   = 0;
 290         }
 291     }
 292     else
 293     {
 294         for (i = 0; i < nfiles; i++)
 295         {
 296             settime[i] = readtime[i];
 297         }
 298     }
 299
 300     if (bSort && (nfiles > 1))
 301     {
 302         sort_files(fnms, settime, nfiles);
 303     }
 304     else
 305     {
 306         fprintf(stderr, "Sorting disabled.\n");
 307     }
 308
 309
 310     /* Write out the new order and start times */
 311     fprintf(stderr, "\nSummary of files and start times used:\n\n"
 312             "          File                Start time\n"
 313             "-----------------------------------------\n");
 314     for (i = 0; i < nfiles; i++)
 315     {
 316         switch (cont_type[i])
 317         {
 318             case TIME_EXPLICIT:
 319                 fprintf(stderr, "%25s   %10.3f\n", fnms[i], settime[i]);
 320                 break;
 321             case TIME_CONTINUE:
 322                 fprintf(stderr, "%25s        Continue from end of last file\n", fnms[i]);
 323                 break;
 324             case TIME_LAST:
 325                 fprintf(stderr, "%25s        Change by same amount as last file\n", fnms[i]);
 326                 break;
 327         }
 328     }
 329     fprintf(stderr, "\n");
 330
 331     settime[nfiles]   = FLT_MAX;
 332     cont_type[nfiles] = TIME_EXPLICIT;
 333     readtime[nfiles]  = FLT_MAX;
 334 }
 335
 336
 337 static void copy_ee(t_energy *src, t_energy *dst, int nre)
 338 {
 339     int i;
 340
 341     for (i = 0; i < nre; i++)
 342     {
 343         dst[i].e    = src[i].e;
 344         dst[i].esum = src[i].esum;
 345         dst[i].eav  = src[i].eav;
 346     }
 347 }
 348
 349 static void update_ee(t_energy *lastee, gmx_large_int_t laststep,
 350                       t_energy *startee, gmx_large_int_t startstep,
 351                       t_energy *ee, int step,
 352                       t_energy *outee, int nre)
 353 {
 354     int    i;
 355     double sigmacorr, nom, denom;
 356     double prestart_esum;
 357     double prestart_sigma;
 358
 359     for (i = 0; i < nre; i++)
 360     {
 361         outee[i].e = ee[i].e;
 362         /* add statistics from earlier file if present */
 363         if (laststep > 0)
 364         {
 365             outee[i].esum = lastee[i].esum+ee[i].esum;
 366             nom           = (lastee[i].esum*(step+1)-ee[i].esum*(laststep));
 367             denom         = ((step+1.0)*(laststep)*(step+1.0+laststep));
 368             sigmacorr     = nom*nom/denom;
 369             outee[i].eav  = lastee[i].eav+ee[i].eav+sigmacorr;
 370         }
 371         else
 372         {
 373             /* otherwise just copy to output */
 374             outee[i].esum = ee[i].esum;
 375             outee[i].eav  = ee[i].eav;
 376         }
 377
 378         /* if we didnt start to write at the first frame
 379          * we must compensate the statistics for this
 380          * there are laststep frames in the earlier file
 381          * and step+1 frames in this one.
 382          */
 383         if (startstep > 0)
 384         {
 385             gmx_large_int_t q = laststep+step;
 386             gmx_large_int_t p = startstep+1;
 387             prestart_esum  = startee[i].esum-startee[i].e;
 388             sigmacorr      = prestart_esum-(p-1)*startee[i].e;
 389             prestart_sigma = startee[i].eav-
 390                 sigmacorr*sigmacorr/(p*(p-1));
 391             sigmacorr = prestart_esum/(p-1)-
 392                 outee[i].esum/(q);
 393             outee[i].esum -= prestart_esum;
 394             if (q-p+1 > 0)
 395             {
 396                 outee[i].eav = outee[i].eav-prestart_sigma-
 397                     sigmacorr*sigmacorr*((p-1)*q)/(q-p+1);
 398             }
 399         }
 400
 401         if ((outee[i].eav/(laststep+step+1)) < (GMX_REAL_EPS))
 402         {
 403             outee[i].eav = 0;
 404         }
 405     }
 406 }
 407
 408 static void update_ee_sum(int nre,
 409                           gmx_large_int_t *ee_sum_step,
 410                           gmx_large_int_t *ee_sum_nsteps,
 411                           gmx_large_int_t *ee_sum_nsum,
 412                           t_energy *ee_sum,
 413                           t_enxframe *fr, int out_step)
 414 {
 415     gmx_large_int_t nsteps, nsum, fr_nsum;
 416     int             i;
 417
 418     nsteps = *ee_sum_nsteps;
 419     nsum   = *ee_sum_nsum;
 420
 421     fr_nsum = fr->nsum;
 422     if (fr_nsum == 0)
 423     {
 424         fr_nsum = 1;
 425     }
 426
 427     if (nsteps == 0)
 428     {
 429         if (fr_nsum == 1)
 430         {
 431             for (i = 0; i < nre; i++)
 432             {
 433                 ee_sum[i].esum = fr->ener[i].e;
 434                 ee_sum[i].eav  = 0;
 435             }
 436         }
 437         else
 438         {
 439             for (i = 0; i < nre; i++)
 440             {
 441                 ee_sum[i].esum = fr->ener[i].esum;
 442                 ee_sum[i].eav  = fr->ener[i].eav;
 443             }
 444         }
 445         nsteps = fr->nsteps;
 446         nsum   = fr_nsum;
 447     }
 448     else if (out_step + *ee_sum_nsum - *ee_sum_step == nsteps + fr->nsteps)
 449     {
 450         if (fr_nsum == 1)
 451         {
 452             for (i = 0; i < nre; i++)
 453             {
 454                 ee_sum[i].eav  +=
 455                     dsqr(ee_sum[i].esum/nsum
 456                          - (ee_sum[i].esum + fr->ener[i].e)/(nsum + 1))*nsum*(nsum + 1);
 457                 ee_sum[i].esum += fr->ener[i].e;
 458             }
 459         }
 460         else
 461         {
 462             for (i = 0; i < fr->nre; i++)
 463             {
 464                 ee_sum[i].eav  +=
 465                     fr->ener[i].eav +
 466                     dsqr(ee_sum[i].esum/nsum
 467                          - (ee_sum[i].esum + fr->ener[i].esum)/(nsum + fr->nsum))*
 468                     nsum*(nsum + fr->nsum)/(double)fr->nsum;
 469                 ee_sum[i].esum += fr->ener[i].esum;
 470             }
 471         }
 472         nsteps += fr->nsteps;
 473         nsum   += fr_nsum;
 474     }
 475     else
 476     {
 477         if (fr->nsum != 0)
 478         {
 479             fprintf(stderr, "\nWARNING: missing energy sums at time %f\n", fr->t);
 480         }
 481         nsteps = 0;
 482         nsum   = 0;
 483     }
 484
 485     *ee_sum_step   = out_step;
 486     *ee_sum_nsteps = nsteps;
 487     *ee_sum_nsum   = nsum;
 488 }
 489
 490 int gmx_eneconv(int argc, char *argv[])
 491 {
 492     const char     *desc[] = {
 493         "With [IT]multiple files[it] specified for the [TT]-f[tt] option:[BR]",
 494         "Concatenates several energy files in sorted order.",
 495         "In the case of double time frames, the one",
 496         "in the later file is used. By specifying [TT]-settime[tt] you will be",
 497         "asked for the start time of each file. The input files are taken",
 498         "from the command line,",
 499         "such that the command [TT]eneconv -f *.edr -o fixed.edr[tt] should do",
 500         "the trick. [PAR]",
 501         "With [IT]one file[it] specified for [TT]-f[tt]:[BR]",
 502         "Reads one energy file and writes another, applying the [TT]-dt[tt],",
 503         "[TT]-offset[tt], [TT]-t0[tt] and [TT]-settime[tt] options and",
 504         "converting to a different format if necessary (indicated by file",
 505         "extentions).[PAR]",
 506         "[TT]-settime[tt] is applied first, then [TT]-dt[tt]/[TT]-offset[tt]",
 507         "followed by [TT]-b[tt] and [TT]-e[tt] to select which frames to write."
 508     };
 509     const char     *bugs[] = {
 510         "When combining trajectories the sigma and E^2 (necessary for statistics) are not updated correctly. Only the actual energy is correct. One thus has to compute statistics in another way."
 511     };
 512     ener_file_t     in  = NULL, out = NULL;
 513     gmx_enxnm_t    *enm = NULL;
 514 #if 0
 515     ener_file_t     in, out = NULL;
 516     gmx_enxnm_t    *enm = NULL;
 517 #endif
 518     t_enxframe     *fr, *fro;
 519     gmx_large_int_t ee_sum_step = 0, ee_sum_nsteps, ee_sum_nsum;
 520     t_energy       *ee_sum;
 521     gmx_large_int_t lastfilestep, laststep, startstep, startstep_file = 0;
 522     int             noutfr;
 523     int             nre, nremax, this_nre, nfile, f, i, j, kkk, nset, *set = NULL;
 524     double          last_t;
 525     char          **fnms;
 526     real           *readtime, *settime, timestep, t1, tadjust;
 527     char            inputstring[STRLEN], *chptr, buf[22], buf2[22], buf3[22];
 528     gmx_bool        ok;
 529     int            *cont_type;
 530     gmx_bool        bNewFile, bFirst, bNewOutput;
 531     output_env_t    oenv;
 532     gmx_bool        warned_about_dh = FALSE;
 533     t_enxblock     *blocks          = NULL;
 534     int             nblocks         = 0;
 535     int             nblocks_alloc   = 0;
 536
 537     t_filenm        fnm[] = {
 538         { efEDR, "-f", NULL,    ffRDMULT },
 539         { efEDR, "-o", "fixed", ffWRITE  },
 540     };
 541
 542 #define NFILE asize(fnm)
 543     gmx_bool         bWrite;
 544     static real      delta_t   = 0.0, toffset = 0, scalefac = 1;
 545     static gmx_bool  bSetTime  = FALSE;
 546     static gmx_bool  bSort     = TRUE, bError = TRUE;
 547     static real      begin     = -1;
 548     static real      end       = -1;
 549     gmx_bool         remove_dh = FALSE;
 550
 551     t_pargs          pa[] = {
 552         { "-b",        FALSE, etREAL, {&begin},
 553           "First time to use"},
 554         { "-e",        FALSE, etREAL, {&end},
 555           "Last time to use"},
 556         { "-dt",       FALSE, etREAL, {&delta_t},
 557           "Only write out frame when t MOD dt = offset" },
 558         { "-offset",   FALSE, etREAL, {&toffset},
 559           "Time offset for [TT]-dt[tt] option" },
 560         { "-settime",  FALSE, etBOOL, {&bSetTime},
 561           "Change starting time interactively" },
 562         { "-sort",     FALSE, etBOOL, {&bSort},
 563           "Sort energy files (not frames)"},
 564         { "-rmdh",     FALSE, etBOOL, {&remove_dh},
 565           "Remove free energy block data" },
 566         { "-scalefac", FALSE, etREAL, {&scalefac},
 567           "Multiply energy component by this factor" },
 568         { "-error",    FALSE, etBOOL, {&bError},
 569           "Stop on errors in the file" }
 570     };
 571
 572     if (!parse_common_args(&argc, argv, PCA_BE_NICE, NFILE, fnm, asize(pa),
 573                            pa, asize(desc), desc, asize(bugs), bugs, &oenv))
 574     {
 575         return 0;
 576     }
 577     tadjust  = 0;
 578     nremax   = 0;
 579     nset     = 0;
 580     timestep = 0.0;
 581     snew(fnms, argc);
 582     nfile        = 0;
 583     lastfilestep = 0;
 584     laststep     = startstep = 0;
 585
 586     nfile = opt2fns(&fnms, "-f", NFILE, fnm);
 587
 588     if (!nfile)
 589     {
 590         gmx_fatal(FARGS, "No input files!");
 591     }
 592
 593     snew(settime, nfile+1);
 594     snew(readtime, nfile+1);
 595     snew(cont_type, nfile+1);
 596
 597     nre = scan_ene_files(fnms, nfile, readtime, &timestep, &nremax);
 598     edit_files(fnms, nfile, readtime, settime, cont_type, bSetTime, bSort);
 599
 600     ee_sum_nsteps = 0;
 601     ee_sum_nsum   = 0;
 602     snew(ee_sum, nremax);
 603
 604     snew(fr, 1);
 605     snew(fro, 1);
 606     fro->t   = -1e20;
 607     fro->nre = nre;
 608     snew(fro->ener, nremax);
 609
 610     noutfr = 0;
 611     bFirst = TRUE;
 612
 613     last_t = fro->t;
 614     for (f = 0; f < nfile; f++)
 615     {
 616         bNewFile   = TRUE;
 617         bNewOutput = TRUE;
 618         in         = open_enx(fnms[f], "r");
 619         enm        = NULL;
 620         do_enxnms(in, &this_nre, &enm);
 621         if (f == 0)
 622         {
 623             if (scalefac != 1)
 624             {
 625                 set = select_it(nre, enm, &nset);
 626             }
 627
 628             /* write names to the output file */
 629             out = open_enx(opt2fn("-o", NFILE, fnm), "w");
 630             do_enxnms(out, &nre, &enm);
 631         }
 632
 633         /* start reading from the next file */
 634         while ((fro->t <= (settime[f+1] + GMX_REAL_EPS)) &&
 635                do_enx(in, fr))
 636         {
 637             if (bNewFile)
 638             {
 639                 startstep_file = fr->step;
 640                 tadjust        = settime[f] - fr->t;
 641                 if (cont_type[f+1] == TIME_LAST)
 642                 {
 643                     settime[f+1]   = readtime[f+1]-readtime[f]+settime[f];
 644                     cont_type[f+1] = TIME_EXPLICIT;
 645                 }
 646                 bNewFile = FALSE;
 647             }
 648
 649             if (tadjust + fr->t <= last_t)
 650             {
 651                 /* Skip this frame, since we already have it / past it */
 652                 if (debug)
 653                 {
 654                     fprintf(debug, "fr->step %s, fr->t %.4f\n",
 655                             gmx_step_str(fr->step, buf), fr->t);
 656                     fprintf(debug, "tadjust %12.6e + fr->t %12.6e <= t %12.6e\n",
 657                             tadjust, fr->t, last_t);
 658                 }
 659                 continue;
 660             }
 661
 662             fro->step = lastfilestep + fr->step - startstep_file;
 663             fro->t    = tadjust  + fr->t;
 664
 665             bWrite = ((begin < 0 || (begin >= 0 && (fro->t >= begin-GMX_REAL_EPS))) &&
 666                       (end  < 0 || (end  >= 0 && (fro->t <= end  +GMX_REAL_EPS))) &&
 667                       (fro->t <= settime[f+1]+0.5*timestep));
 668
 669             if (debug)
 670             {
 671                 fprintf(debug,
 672                         "fr->step %s, fr->t %.4f, fro->step %s fro->t %.4f, w %d\n",
 673                         gmx_step_str(fr->step, buf), fr->t,
 674                         gmx_step_str(fro->step, buf2), fro->t, bWrite);
 675             }
 676
 677             if (bError)
 678             {
 679                 if ((end > 0) && (fro->t > end+GMX_REAL_EPS))
 680                 {
 681                     f = nfile;
 682                     break;
 683                 }
 684             }
 685
 686             if (fro->t >= begin-GMX_REAL_EPS)
 687             {
 688                 if (bFirst)
 689                 {
 690                     bFirst    = FALSE;
 691                     startstep = fr->step;
 692                 }
 693                 if (bWrite)
 694                 {
 695                     update_ee_sum(nre, &ee_sum_step, &ee_sum_nsteps, &ee_sum_nsum, ee_sum,
 696                                   fr, fro->step);
 697                 }
 698             }
 699
 700             /* determine if we should write it */
 701             if (bWrite && (delta_t == 0 || bRmod(fro->t, toffset, delta_t)))
 702             {
 703                 laststep = fro->step;
 704                 last_t   = fro->t;
 705                 if (bNewOutput)
 706                 {
 707                     bNewOutput = FALSE;
 708                     fprintf(stderr, "\nContinue writing frames from t=%g, step=%s\n",
 709                             fro->t, gmx_step_str(fro->step, buf));
 710                 }
 711
 712                 /* Copy the energies */
 713                 for (i = 0; i < nre; i++)
 714                 {
 715                     fro->ener[i].e = fr->ener[i].e;
 716                 }
 717
 718                 fro->nsteps = ee_sum_nsteps;
 719                 fro->dt     = fr->dt;
 720
 721                 if (ee_sum_nsum <= 1)
 722                 {
 723                     fro->nsum = 0;
 724                 }
 725                 else
 726                 {
 727                     fro->nsum = gmx_large_int_to_int(ee_sum_nsum,
 728                                                      "energy average summation");
 729                     /* Copy the energy sums */
 730                     for (i = 0; i < nre; i++)
 731                     {
 732                         fro->ener[i].esum = ee_sum[i].esum;
 733                         fro->ener[i].eav  = ee_sum[i].eav;
 734                     }
 735                 }
 736                 /* We wrote the energies, so reset the counts */
 737                 ee_sum_nsteps = 0;
 738                 ee_sum_nsum   = 0;
 739
 740                 if (scalefac != 1)
 741                 {
 742                     for (kkk = 0; kkk < nset; kkk++)
 743                     {
 744                         fro->ener[set[kkk]].e    *= scalefac;
 745                         if (fro->nsum > 0)
 746                         {
 747                             fro->ener[set[kkk]].eav  *= scalefac*scalefac;
 748                             fro->ener[set[kkk]].esum *= scalefac;
 749                         }
 750                     }
 751                 }
 752                 /* Copy restraint stuff */
 753                 /*fro->ndisre       = fr->ndisre;
 754                    fro->disre_rm3tav = fr->disre_rm3tav;
 755                    fro->disre_rt     = fr->disre_rt;*/
 756                 fro->nblock       = fr->nblock;
 757                 /*fro->nr           = fr->nr;*/
 758                 fro->block        = fr->block;
 759
 760                 /* check if we have blocks with delta_h data and are throwing
 761                    away data */
 762                 if (fro->nblock > 0)
 763                 {
 764                     if (remove_dh)
 765                     {
 766                         int i;
 767                         if (!blocks || nblocks_alloc < fr->nblock)
 768                         {
 769                             /* we pre-allocate the blocks */
 770                             nblocks_alloc = fr->nblock;
 771                             snew(blocks, nblocks_alloc);
 772                         }
 773                         nblocks = 0; /* number of blocks so far */
 774
 775                         for (i = 0; i < fr->nblock; i++)
 776                         {
 777                             if ( (fr->block[i].id != enxDHCOLL) &&
 778                                  (fr->block[i].id != enxDH) &&
 779                                  (fr->block[i].id != enxDHHIST) )
 780                             {
 781                                 /* copy everything verbatim */
 782                                 blocks[nblocks] = fr->block[i];
 783                                 nblocks++;
 784                             }
 785                         }
 786                         /* now set the block pointer to the new blocks */
 787                         fro->nblock = nblocks;
 788                         fro->block  = blocks;
 789                     }
 790                     else if (delta_t > 0)
 791                     {
 792                         if (!warned_about_dh)
 793                         {
 794                             for (i = 0; i < fr->nblock; i++)
 795                             {
 796                                 if (fr->block[i].id == enxDH ||
 797                                     fr->block[i].id == enxDHHIST)
 798                                 {
 799                                     int size;
 800                                     if (fr->block[i].id == enxDH)
 801                                     {
 802                                         size = fr->block[i].sub[2].nr;
 803                                     }
 804                                     else
 805                                     {
 806                                         size = fr->nsteps;
 807                                     }
 808                                     if (size > 0)
 809                                     {
 810                                         printf("\nWARNING: %s contains delta H blocks or histograms for which\n"
 811                                                "         some data is thrown away on a block-by-block basis, where each block\n"
 812                                                "         contains up to %d samples.\n"
 813                                                "         This is almost certainly not what you want.\n"
 814                                                "         Use the -rmdh option to throw all delta H samples away.\n"
 815                                                "         Use g_energy -odh option to extract these samples.\n",
 816                                                fnms[f], size);
 817                                         warned_about_dh = TRUE;
 818                                         break;
 819                                     }
 820                                 }
 821                             }
 822                         }
 823                     }
 824                 }
 825
 826                 do_enx(out, fro);
 827                 if (noutfr % 1000 == 0)
 828                 {
 829                     fprintf(stderr, "Writing frame time %g    ", fro->t);
 830                 }
 831                 noutfr++;
 832             }
 833         }
 834         if (f == nfile)
 835         {
 836             f--;
 837         }
 838         printf("\nLast step written from %s: t %g, step %s\n",
 839                fnms[f], last_t, gmx_step_str(laststep, buf));
 840         lastfilestep = laststep;
 841
 842         /* set the next time from the last in previous file */
 843         if (cont_type[f+1] == TIME_CONTINUE)
 844         {
 845             settime[f+1] = fro->t;
 846             /* in this case we have already written the last frame of
 847              * previous file, so update begin to avoid doubling it
 848              * with the start of the next file
 849              */
 850             begin = fro->t+0.5*timestep;
 851             /* cont_type[f+1]==TIME_EXPLICIT; */
 852         }
 853
 854         if ((fro->t < end) && (f < nfile-1) &&
 855             (fro->t < settime[f+1]-1.5*timestep))
 856         {
 857             fprintf(stderr,
 858                     "\nWARNING: There might be a gap around t=%g\n", fro->t);
 859         }
 860
 861         /* move energies to lastee */
 862         close_enx(in);
 863         free_enxnms(this_nre, enm);
 864
 865         fprintf(stderr, "\n");
 866     }
 867     if (noutfr == 0)
 868     {
 869         fprintf(stderr, "No frames written.\n");
 870     }
 871     else
 872     {
 873         fprintf(stderr, "Last frame written was at step %s, time %f\n",
 874                 gmx_step_str(fro->step, buf), fro->t);
 875         fprintf(stderr, "Wrote %d frames\n", noutfr);
 876     }
 877
 878     return 0;
 879 }