src/gromacs/gmxana/gmx_eneconv.cpp

   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 <cmath>
  40 #include <cstdlib>
  41 #include <cstring>
  42
  43 #include <algorithm>
  44
  45 #include "gromacs/commandline/pargs.h"
  46 #include "gromacs/fileio/enxio.h"
  47 #include "gromacs/fileio/trxio.h"
  48 #include "gromacs/gmxana/gmx_ana.h"
  49 #include "gromacs/legacyheaders/disre.h"
  50 #include "gromacs/legacyheaders/macros.h"
  51 #include "gromacs/legacyheaders/names.h"
  52 #include "gromacs/legacyheaders/typedefs.h"
  53 #include "gromacs/math/vec.h"
  54 #include "gromacs/utility/cstringutil.h"
  55 #include "gromacs/utility/fatalerror.h"
  56 #include "gromacs/utility/smalloc.h"
  57
  58 #define TIME_EXPLICIT 0
  59 #define TIME_CONTINUE 1
  60 #define TIME_LAST     2
  61 #ifndef FLT_MAX
  62 #define FLT_MAX 1e36
  63 #endif
  64
  65 static int *select_it(int nre, gmx_enxnm_t *nm, int *nset)
  66 {
  67     gmx_bool *bE;
  68     int       n, k, j, i;
  69     int      *set;
  70     gmx_bool  bVerbose = TRUE;
  71
  72     if ((getenv("GMX_ENER_VERBOSE")) != NULL)
  73     {
  74         bVerbose = FALSE;
  75     }
  76
  77     fprintf(stderr, "Select the terms you want to scale from the following list\n");
  78     fprintf(stderr, "End your selection with 0\n");
  79
  80     if (bVerbose)
  81     {
  82         for (k = 0; (k < nre); )
  83         {
  84             for (j = 0; (j < 4) && (k < nre); j++, k++)
  85             {
  86                 fprintf(stderr, " %3d=%14s", k+1, nm[k].name);
  87             }
  88             fprintf(stderr, "\n");
  89         }
  90     }
  91
  92     snew(bE, nre);
  93     do
  94     {
  95         if (1 != scanf("%d", &n))
  96         {
  97             gmx_fatal(FARGS, "Cannot read energy term");
  98         }
  99         if ((n > 0) && (n <= nre))
 100         {
 101             bE[n-1] = TRUE;
 102         }
 103     }
 104     while (n != 0);
 105
 106     snew(set, nre);
 107     for (i = (*nset) = 0; (i < nre); i++)
 108     {
 109         if (bE[i])
 110         {
 111             set[(*nset)++] = i;
 112         }
 113     }
 114
 115     sfree(bE);
 116
 117     return set;
 118 }
 119
 120 static void sort_files(char **fnms, real *settime, int nfile)
 121 {
 122     int   i, j, minidx;
 123     real  timeswap;
 124     char *chptr;
 125
 126     for (i = 0; i < nfile; i++)
 127     {
 128         minidx = i;
 129         for (j = i+1; j < nfile; j++)
 130         {
 131             if (settime[j] < settime[minidx])
 132             {
 133                 minidx = j;
 134             }
 135         }
 136         if (minidx != i)
 137         {
 138             timeswap        = settime[i];
 139             settime[i]      = settime[minidx];
 140             settime[minidx] = timeswap;
 141             chptr           = fnms[i];
 142             fnms[i]         = fnms[minidx];
 143             fnms[minidx]    = chptr;
 144         }
 145     }
 146 }
 147
 148
 149 static int scan_ene_files(char **fnms, int nfiles,
 150                           real *readtime, real *timestep, int *nremax)
 151 {
 152     /* Check number of energy terms and start time of all files */
 153     int          f, nre, nremin = 0, nresav = 0;
 154     ener_file_t  in;
 155     real         t1, t2;
 156     char         inputstring[STRLEN];
 157     gmx_enxnm_t *enm;
 158     t_enxframe  *fr;
 159
 160     snew(fr, 1);
 161
 162     for (f = 0; f < nfiles; f++)
 163     {
 164         in  = open_enx(fnms[f], "r");
 165         enm = NULL;
 166         do_enxnms(in, &nre, &enm);
 167
 168         if (f == 0)
 169         {
 170             nresav  = nre;
 171             nremin  = nre;
 172             *nremax = nre;
 173             do_enx(in, fr);
 174             t1 = fr->t;
 175             do_enx(in, fr);
 176             t2          = fr->t;
 177             *timestep   = t2-t1;
 178             readtime[f] = t1;
 179             close_enx(in);
 180         }
 181         else
 182         {
 183             nremin  = std::min(nremin, fr->nre);
 184             *nremax = std::max(*nremax, fr->nre);
 185             if (nre != nresav)
 186             {
 187                 fprintf(stderr,
 188                         "Energy files don't match, different number of energies:\n"
 189                         " %s: %d\n %s: %d\n", fnms[f-1], nresav, fnms[f], fr->nre);
 190                 fprintf(stderr,
 191                         "\nContinue conversion using only the first %d terms (n/y)?\n"
 192                         "(you should be sure that the energy terms match)\n", nremin);
 193                 if (NULL == fgets(inputstring, STRLEN-1, stdin))
 194                 {
 195                     gmx_fatal(FARGS, "Error reading user input");
 196                 }
 197                 if (inputstring[0] != 'y' && inputstring[0] != 'Y')
 198                 {
 199                     fprintf(stderr, "Will not convert\n");
 200                     exit(0);
 201                 }
 202                 nresav = fr->nre;
 203             }
 204             do_enx(in, fr);
 205             readtime[f] = fr->t;
 206             close_enx(in);
 207         }
 208         fprintf(stderr, "\n");
 209         free_enxnms(nre, enm);
 210     }
 211
 212     free_enxframe(fr);
 213     sfree(fr);
 214
 215     return nremin;
 216 }
 217
 218
 219 static void edit_files(char **fnms, int nfiles, real *readtime,
 220                        real *settime, int *cont_type, gmx_bool bSetTime, gmx_bool bSort)
 221 {
 222     int      i;
 223     gmx_bool ok;
 224     char     inputstring[STRLEN], *chptr;
 225
 226     if (bSetTime)
 227     {
 228         if (nfiles == 1)
 229         {
 230             fprintf(stderr, "\n\nEnter the new start time:\n\n");
 231         }
 232         else
 233         {
 234             fprintf(stderr, "\n\nEnter the new start time for each file.\n"
 235                     "There are two special options, both disables sorting:\n\n"
 236                     "c (continue) - The start time is taken from the end\n"
 237                     "of the previous file. Use it when your continuation run\n"
 238                     "restarts with t=0 and there is no overlap.\n\n"
 239                     "l (last) - The time in this file will be changed the\n"
 240                     "same amount as in the previous. Use it when the time in the\n"
 241                     "new run continues from the end of the previous one,\n"
 242                     "since this takes possible overlap into account.\n\n");
 243         }
 244
 245         fprintf(stderr, "          File             Current start       New start\n"
 246                 "---------------------------------------------------------\n");
 247
 248         for (i = 0; i < nfiles; i++)
 249         {
 250             fprintf(stderr, "%25s   %10.3f             ", fnms[i], readtime[i]);
 251             ok = FALSE;
 252             do
 253             {
 254                 if (NULL == fgets(inputstring, STRLEN-1, stdin))
 255                 {
 256                     gmx_fatal(FARGS, "Error reading user input");
 257                 }
 258                 inputstring[std::strlen(inputstring)-1] = 0;
 259
 260                 if (inputstring[0] == 'c' || inputstring[0] == 'C')
 261                 {
 262                     cont_type[i] = TIME_CONTINUE;
 263                     bSort        = FALSE;
 264                     ok           = TRUE;
 265                     settime[i]   = FLT_MAX;
 266                 }
 267                 else if (inputstring[0] == 'l' ||
 268                          inputstring[0] == 'L')
 269                 {
 270                     cont_type[i] = TIME_LAST;
 271                     bSort        = FALSE;
 272                     ok           = TRUE;
 273                     settime[i]   = FLT_MAX;
 274                 }
 275                 else
 276                 {
 277                     settime[i] = strtod(inputstring, &chptr);
 278                     if (chptr == inputstring)
 279                     {
 280                         fprintf(stderr, "Try that again: ");
 281                     }
 282                     else
 283                     {
 284                         cont_type[i] = TIME_EXPLICIT;
 285                         ok           = TRUE;
 286                     }
 287                 }
 288             }
 289             while (!ok);
 290         }
 291         if (cont_type[0] != TIME_EXPLICIT)
 292         {
 293             cont_type[0] = TIME_EXPLICIT;
 294             settime[0]   = 0;
 295         }
 296     }
 297     else
 298     {
 299         for (i = 0; i < nfiles; i++)
 300         {
 301             settime[i] = readtime[i];
 302         }
 303     }
 304
 305     if (bSort && (nfiles > 1))
 306     {
 307         sort_files(fnms, settime, nfiles);
 308     }
 309     else
 310     {
 311         fprintf(stderr, "Sorting disabled.\n");
 312     }
 313
 314
 315     /* Write out the new order and start times */
 316     fprintf(stderr, "\nSummary of files and start times used:\n\n"
 317             "          File                Start time\n"
 318             "-----------------------------------------\n");
 319     for (i = 0; i < nfiles; i++)
 320     {
 321         switch (cont_type[i])
 322         {
 323             case TIME_EXPLICIT:
 324                 fprintf(stderr, "%25s   %10.3f\n", fnms[i], settime[i]);
 325                 break;
 326             case TIME_CONTINUE:
 327                 fprintf(stderr, "%25s        Continue from end of last file\n", fnms[i]);
 328                 break;
 329             case TIME_LAST:
 330                 fprintf(stderr, "%25s        Change by same amount as last file\n", fnms[i]);
 331                 break;
 332         }
 333     }
 334     fprintf(stderr, "\n");
 335
 336     settime[nfiles]   = FLT_MAX;
 337     cont_type[nfiles] = TIME_EXPLICIT;
 338     readtime[nfiles]  = FLT_MAX;
 339 }
 340
 341
 342 static void copy_ee(t_energy *src, t_energy *dst, int nre)
 343 {
 344     int i;
 345
 346     for (i = 0; i < nre; i++)
 347     {
 348         dst[i].e    = src[i].e;
 349         dst[i].esum = src[i].esum;
 350         dst[i].eav  = src[i].eav;
 351     }
 352 }
 353
 354 static void update_ee(t_energy *lastee, gmx_int64_t laststep,
 355                       t_energy *startee, gmx_int64_t startstep,
 356                       t_energy *ee, int step,
 357                       t_energy *outee, int nre)
 358 {
 359     int    i;
 360     double sigmacorr, nom, denom;
 361     double prestart_esum;
 362     double prestart_sigma;
 363
 364     for (i = 0; i < nre; i++)
 365     {
 366         outee[i].e = ee[i].e;
 367         /* add statistics from earlier file if present */
 368         if (laststep > 0)
 369         {
 370             outee[i].esum = lastee[i].esum+ee[i].esum;
 371             nom           = (lastee[i].esum*(step+1)-ee[i].esum*(laststep));
 372             denom         = ((step+1.0)*(laststep)*(step+1.0+laststep));
 373             sigmacorr     = nom*nom/denom;
 374             outee[i].eav  = lastee[i].eav+ee[i].eav+sigmacorr;
 375         }
 376         else
 377         {
 378             /* otherwise just copy to output */
 379             outee[i].esum = ee[i].esum;
 380             outee[i].eav  = ee[i].eav;
 381         }
 382
 383         /* if we didnt start to write at the first frame
 384          * we must compensate the statistics for this
 385          * there are laststep frames in the earlier file
 386          * and step+1 frames in this one.
 387          */
 388         if (startstep > 0)
 389         {
 390             gmx_int64_t q = laststep+step;
 391             gmx_int64_t p = startstep+1;
 392             prestart_esum  = startee[i].esum-startee[i].e;
 393             sigmacorr      = prestart_esum-(p-1)*startee[i].e;
 394             prestart_sigma = startee[i].eav-
 395                 sigmacorr*sigmacorr/(p*(p-1));
 396             sigmacorr = prestart_esum/(p-1)-
 397                 outee[i].esum/(q);
 398             outee[i].esum -= prestart_esum;
 399             if (q-p+1 > 0)
 400             {
 401                 outee[i].eav = outee[i].eav-prestart_sigma-
 402                     sigmacorr*sigmacorr*((p-1)*q)/(q-p+1);
 403             }
 404         }
 405
 406         if ((outee[i].eav/(laststep+step+1)) < (GMX_REAL_EPS))
 407         {
 408             outee[i].eav = 0;
 409         }
 410     }
 411 }
 412
 413 static void update_ee_sum(int nre,
 414                           gmx_int64_t *ee_sum_step,
 415                           gmx_int64_t *ee_sum_nsteps,
 416                           gmx_int64_t *ee_sum_nsum,
 417                           t_energy *ee_sum,
 418                           t_enxframe *fr, int out_step)
 419 {
 420     gmx_int64_t     nsteps, nsum, fr_nsum;
 421     int             i;
 422
 423     nsteps = *ee_sum_nsteps;
 424     nsum   = *ee_sum_nsum;
 425
 426     fr_nsum = fr->nsum;
 427     if (fr_nsum == 0)
 428     {
 429         fr_nsum = 1;
 430     }
 431
 432     if (nsteps == 0)
 433     {
 434         if (fr_nsum == 1)
 435         {
 436             for (i = 0; i < nre; i++)
 437             {
 438                 ee_sum[i].esum = fr->ener[i].e;
 439                 ee_sum[i].eav  = 0;
 440             }
 441         }
 442         else
 443         {
 444             for (i = 0; i < nre; i++)
 445             {
 446                 ee_sum[i].esum = fr->ener[i].esum;
 447                 ee_sum[i].eav  = fr->ener[i].eav;
 448             }
 449         }
 450         nsteps = fr->nsteps;
 451         nsum   = fr_nsum;
 452     }
 453     else if (out_step + *ee_sum_nsum - *ee_sum_step == nsteps + fr->nsteps)
 454     {
 455         if (fr_nsum == 1)
 456         {
 457             for (i = 0; i < nre; i++)
 458             {
 459                 ee_sum[i].eav  +=
 460                     dsqr(ee_sum[i].esum/nsum
 461                          - (ee_sum[i].esum + fr->ener[i].e)/(nsum + 1))*nsum*(nsum + 1);
 462                 ee_sum[i].esum += fr->ener[i].e;
 463             }
 464         }
 465         else
 466         {
 467             for (i = 0; i < fr->nre; i++)
 468             {
 469                 ee_sum[i].eav  +=
 470                     fr->ener[i].eav +
 471                     dsqr(ee_sum[i].esum/nsum
 472                          - (ee_sum[i].esum + fr->ener[i].esum)/(nsum + fr->nsum))*
 473                     nsum*(nsum + fr->nsum)/static_cast<double>(fr->nsum);
 474                 ee_sum[i].esum += fr->ener[i].esum;
 475             }
 476         }
 477         nsteps += fr->nsteps;
 478         nsum   += fr_nsum;
 479     }
 480     else
 481     {
 482         if (fr->nsum != 0)
 483         {
 484             fprintf(stderr, "\nWARNING: missing energy sums at time %f\n", fr->t);
 485         }
 486         nsteps = 0;
 487         nsum   = 0;
 488     }
 489
 490     *ee_sum_step   = out_step;
 491     *ee_sum_nsteps = nsteps;
 492     *ee_sum_nsum   = nsum;
 493 }
 494
 495 int gmx_eneconv(int argc, char *argv[])
 496 {
 497     const char     *desc[] = {
 498         "With [IT]multiple files[it] specified for the [TT]-f[tt] option:[PAR]",
 499         "Concatenates several energy files in sorted order.",
 500         "In the case of double time frames, the one",
 501         "in the later file is used. By specifying [TT]-settime[tt] you will be",
 502         "asked for the start time of each file. The input files are taken",
 503         "from the command line,",
 504         "such that the command [TT]gmx eneconv -f *.edr -o fixed.edr[tt] should do",
 505         "the trick. [PAR]",
 506         "With [IT]one file[it] specified for [TT]-f[tt]:[PAR]",
 507         "Reads one energy file and writes another, applying the [TT]-dt[tt],",
 508         "[TT]-offset[tt], [TT]-t0[tt] and [TT]-settime[tt] options and",
 509         "converting to a different format if necessary (indicated by file",
 510         "extentions).[PAR]",
 511         "[TT]-settime[tt] is applied first, then [TT]-dt[tt]/[TT]-offset[tt]",
 512         "followed by [TT]-b[tt] and [TT]-e[tt] to select which frames to write."
 513     };
 514     const char     *bugs[] = {
 515         "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."
 516     };
 517     ener_file_t     in  = NULL, out = NULL;
 518     gmx_enxnm_t    *enm = NULL;
 519 #if 0
 520     ener_file_t     in, out = NULL;
 521     gmx_enxnm_t    *enm = NULL;
 522 #endif
 523     t_enxframe     *fr, *fro;
 524     gmx_int64_t     ee_sum_step = 0, ee_sum_nsteps, ee_sum_nsum;
 525     t_energy       *ee_sum;
 526     gmx_int64_t     lastfilestep, laststep, startstep_file = 0;
 527     int             noutfr;
 528     int             nre, nremax, this_nre, nfile, f, i, kkk, nset, *set = NULL;
 529     double          last_t;
 530     char          **fnms;
 531     real           *readtime, *settime, timestep, tadjust;
 532     char            buf[22], buf2[22];
 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     = 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                 }
 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 }