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