src/gromacs/gmxana/gmx_eneconv.c

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