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