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,2016,2017,2018, 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/listed-forces/disre.h"
  50 #include "gromacs/math/functions.h"
  51 #include "gromacs/math/vec.h"
  52 #include "gromacs/mdtypes/md_enums.h"
  53 #include "gromacs/trajectory/energyframe.h"
  54 #include "gromacs/utility/arrayref.h"
  55 #include "gromacs/utility/arraysize.h"
  56 #include "gromacs/utility/cstringutil.h"
  57 #include "gromacs/utility/fatalerror.h"
  58 #include "gromacs/utility/int64_to_int.h"
  59 #include "gromacs/utility/smalloc.h"
  60 #include "gromacs/utility/strconvert.h"
  61
  62 #define TIME_EXPLICIT 0
  63 #define TIME_CONTINUE 1
  64 #define TIME_LAST     2
  65 #ifndef FLT_MAX
  66 #define FLT_MAX 1e36
  67 #endif
  68
  69 static int *select_it(int nre, gmx_enxnm_t *nm, int *nset)
  70 {
  71     gmx_bool *bE;
  72     int       n, k, j, i;
  73     int      *set;
  74     gmx_bool  bVerbose = TRUE;
  75
  76     if ((getenv("GMX_ENER_VERBOSE")) != nullptr)
  77     {
  78         bVerbose = FALSE;
  79     }
  80
  81     fprintf(stderr, "Select the terms you want to scale from the following list\n");
  82     fprintf(stderr, "End your selection with 0\n");
  83
  84     if (bVerbose)
  85     {
  86         for (k = 0; (k < nre); )
  87         {
  88             for (j = 0; (j < 4) && (k < nre); j++, k++)
  89             {
  90                 fprintf(stderr, " %3d=%14s", k+1, nm[k].name);
  91             }
  92             fprintf(stderr, "\n");
  93         }
  94     }
  95
  96     snew(bE, nre);
  97     do
  98     {
  99         if (1 != scanf("%d", &n))
 100         {
 101             gmx_fatal(FARGS, "Cannot read energy term");
 102         }
 103         if ((n > 0) && (n <= nre))
 104         {
 105             bE[n-1] = TRUE;
 106         }
 107     }
 108     while (n != 0);
 109
 110     snew(set, nre);
 111     for (i = (*nset) = 0; (i < nre); i++)
 112     {
 113         if (bE[i])
 114         {
 115             set[(*nset)++] = i;
 116         }
 117     }
 118
 119     sfree(bE);
 120
 121     return set;
 122 }
 123
 124 static void sort_files(gmx::ArrayRef<std::string> files, real *settime)
 125 {
 126     for (gmx::index i = 0; i < files.size(); i++)
 127     {
 128         gmx::index minidx = i;
 129         for (gmx::index j = i + 1; j < files.size(); j++)
 130         {
 131             if (settime[j] < settime[minidx])
 132             {
 133                 minidx = j;
 134             }
 135         }
 136         if (minidx != i)
 137         {
 138             real timeswap   = settime[i];
 139             settime[i]      = settime[minidx];
 140             settime[minidx] = timeswap;
 141             std::string tmp = files[i];
 142             files[i]        = files[minidx];
 143             files[minidx]   = tmp;
 144         }
 145     }
 146 }
 147
 148
 149 static int scan_ene_files(const std::vector<std::string> &files,
 150                           real *readtime, real *timestep, int *nremax)
 151 {
 152     /* Check number of energy terms and start time of all files */
 153     int          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 (size_t f = 0; f < files.size(); f++)
 163     {
 164         in  = open_enx(files[f].c_str(), "r");
 165         enm = nullptr;
 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", files[f - 1].c_str(), nresav, files[f].c_str(), 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 (nullptr == 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(gmx::ArrayRef<std::string> files, real *readtime,
 220                        real *settime, int *cont_type, gmx_bool bSetTime, gmx_bool bSort)
 221 {
 222     gmx_bool ok;
 223     char     inputstring[STRLEN], *chptr;
 224
 225     if (bSetTime)
 226     {
 227         if (files.size() == 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 (gmx::index i = 0; i < files.size(); i++)
 248         {
 249             fprintf(stderr, "%25s   %10.3f             ", files[i].c_str(), readtime[i]);
 250             ok = FALSE;
 251             do
 252             {
 253                 if (nullptr == fgets(inputstring, STRLEN-1, stdin))
 254                 {
 255                     gmx_fatal(FARGS, "Error reading user input");
 256                 }
 257                 inputstring[std::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 (gmx::index i = 0; i < files.size(); i++)
 299         {
 300             settime[i] = readtime[i];
 301         }
 302     }
 303
 304     if (bSort && files.size() > 1)
 305     {
 306         sort_files(files, settime);
 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 (gmx::index i = 0; i < files.size(); i++)
 319     {
 320         switch (cont_type[i])
 321         {
 322             case TIME_EXPLICIT:
 323                 fprintf(stderr, "%25s   %10.3f\n", files[i].c_str(), settime[i]);
 324                 break;
 325             case TIME_CONTINUE:
 326                 fprintf(stderr, "%25s        Continue from end of last file\n", files[i].c_str());
 327                 break;
 328             case TIME_LAST:
 329                 fprintf(stderr, "%25s        Change by same amount as last file\n", files[i].c_str());
 330                 break;
 331         }
 332     }
 333     fprintf(stderr, "\n");
 334
 335     settime[files.size()]   = FLT_MAX;
 336     cont_type[files.size()] = TIME_EXPLICIT;
 337     readtime[files.size()]  = FLT_MAX;
 338 }
 339
 340
 341 static void update_ee_sum(int nre,
 342                           int64_t *ee_sum_step,
 343                           int64_t *ee_sum_nsteps,
 344                           int64_t *ee_sum_nsum,
 345                           t_energy *ee_sum,
 346                           t_enxframe *fr, int out_step)
 347 {
 348     int64_t         nsteps, nsum, fr_nsum;
 349     int             i;
 350
 351     nsteps = *ee_sum_nsteps;
 352     nsum   = *ee_sum_nsum;
 353
 354     fr_nsum = fr->nsum;
 355     if (fr_nsum == 0)
 356     {
 357         fr_nsum = 1;
 358     }
 359
 360     if (nsteps == 0)
 361     {
 362         if (fr_nsum == 1)
 363         {
 364             for (i = 0; i < nre; i++)
 365             {
 366                 ee_sum[i].esum = fr->ener[i].e;
 367                 ee_sum[i].eav  = 0;
 368             }
 369         }
 370         else
 371         {
 372             for (i = 0; i < nre; i++)
 373             {
 374                 ee_sum[i].esum = fr->ener[i].esum;
 375                 ee_sum[i].eav  = fr->ener[i].eav;
 376             }
 377         }
 378         nsteps = fr->nsteps;
 379         nsum   = fr_nsum;
 380     }
 381     else if (out_step + *ee_sum_nsum - *ee_sum_step == nsteps + fr->nsteps)
 382     {
 383         if (fr_nsum == 1)
 384         {
 385             for (i = 0; i < nre; i++)
 386             {
 387                 ee_sum[i].eav  +=
 388                     gmx::square(ee_sum[i].esum/nsum
 389                                 - (ee_sum[i].esum + fr->ener[i].e)/(nsum + 1))*nsum*(nsum + 1);
 390                 ee_sum[i].esum += fr->ener[i].e;
 391             }
 392         }
 393         else
 394         {
 395             for (i = 0; i < fr->nre; i++)
 396             {
 397                 ee_sum[i].eav  +=
 398                     fr->ener[i].eav +
 399                     gmx::square(ee_sum[i].esum/nsum
 400                                 - (ee_sum[i].esum + fr->ener[i].esum)/(nsum + fr->nsum))*
 401                     nsum*(nsum + fr->nsum)/static_cast<double>(fr->nsum);
 402                 ee_sum[i].esum += fr->ener[i].esum;
 403             }
 404         }
 405         nsteps += fr->nsteps;
 406         nsum   += fr_nsum;
 407     }
 408     else
 409     {
 410         if (fr->nsum != 0)
 411         {
 412             fprintf(stderr, "\nWARNING: missing energy sums at time %f\n", fr->t);
 413         }
 414         nsteps = 0;
 415         nsum   = 0;
 416     }
 417
 418     *ee_sum_step   = out_step;
 419     *ee_sum_nsteps = nsteps;
 420     *ee_sum_nsum   = nsum;
 421 }
 422
 423 int gmx_eneconv(int argc, char *argv[])
 424 {
 425     const char       *desc[] = {
 426         "With [IT]multiple files[it] specified for the [TT]-f[tt] option:[PAR]",
 427         "Concatenates several energy files in sorted order.",
 428         "In the case of double time frames, the one",
 429         "in the later file is used. By specifying [TT]-settime[tt] you will be",
 430         "asked for the start time of each file. The input files are taken",
 431         "from the command line,",
 432         "such that the command [TT]gmx eneconv -f *.edr -o fixed.edr[tt] should do",
 433         "the trick. [PAR]",
 434         "With [IT]one file[it] specified for [TT]-f[tt]:[PAR]",
 435         "Reads one energy file and writes another, applying the [TT]-dt[tt],",
 436         "[TT]-offset[tt], [TT]-t0[tt] and [TT]-settime[tt] options and",
 437         "converting to a different format if necessary (indicated by file",
 438         "extentions).[PAR]",
 439         "[TT]-settime[tt] is applied first, then [TT]-dt[tt]/[TT]-offset[tt]",
 440         "followed by [TT]-b[tt] and [TT]-e[tt] to select which frames to write."
 441     };
 442     const char       *bugs[] = {
 443         "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."
 444     };
 445     ener_file_t       in  = nullptr, out = nullptr;
 446     gmx_enxnm_t      *enm = nullptr;
 447 #if 0
 448     ener_file_t       in, out = NULL;
 449     gmx_enxnm_t      *enm = NULL;
 450 #endif
 451     t_enxframe       *fr, *fro;
 452     int64_t           ee_sum_step = 0, ee_sum_nsteps, ee_sum_nsum;
 453     t_energy         *ee_sum;
 454     int64_t           lastfilestep, laststep, startstep_file = 0;
 455     int               noutfr;
 456     int               nre, nremax, this_nre, i, kkk, nset, *set = nullptr;
 457     double            last_t;
 458     real             *readtime, *settime, timestep, tadjust;
 459     char              buf[22], buf2[22];
 460     int              *cont_type;
 461     gmx_bool          bNewFile, bFirst, bNewOutput;
 462     gmx_output_env_t *oenv;
 463     gmx_bool          warned_about_dh = FALSE;
 464     t_enxblock       *blocks          = nullptr;
 465     int               nblocks         = 0;
 466     int               nblocks_alloc   = 0;
 467
 468     t_filenm          fnm[] = {
 469         { efEDR, "-f", nullptr,    ffRDMULT },
 470         { efEDR, "-o", "fixed", ffWRITE  },
 471     };
 472
 473 #define NFILE asize(fnm)
 474     gmx_bool         bWrite;
 475     static real      delta_t   = 0.0, toffset = 0, scalefac = 1;
 476     static gmx_bool  bSetTime  = FALSE;
 477     static gmx_bool  bSort     = TRUE, bError = TRUE;
 478     static real      begin     = -1;
 479     static real      end       = -1;
 480     gmx_bool         remove_dh = FALSE;
 481
 482     t_pargs          pa[] = {
 483         { "-b",        FALSE, etREAL, {&begin},
 484           "First time to use"},
 485         { "-e",        FALSE, etREAL, {&end},
 486           "Last time to use"},
 487         { "-dt",       FALSE, etREAL, {&delta_t},
 488           "Only write out frame when t MOD dt = offset" },
 489         { "-offset",   FALSE, etREAL, {&toffset},
 490           "Time offset for [TT]-dt[tt] option" },
 491         { "-settime",  FALSE, etBOOL, {&bSetTime},
 492           "Change starting time interactively" },
 493         { "-sort",     FALSE, etBOOL, {&bSort},
 494           "Sort energy files (not frames)"},
 495         { "-rmdh",     FALSE, etBOOL, {&remove_dh},
 496           "Remove free energy block data" },
 497         { "-scalefac", FALSE, etREAL, {&scalefac},
 498           "Multiply energy component by this factor" },
 499         { "-error",    FALSE, etBOOL, {&bError},
 500           "Stop on errors in the file" }
 501     };
 502
 503     if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa),
 504                            pa, asize(desc), desc, asize(bugs), bugs, &oenv))
 505     {
 506         return 0;
 507     }
 508     tadjust      = 0;
 509     nremax       = 0;
 510     nset         = 0;
 511     timestep     = 0.0;
 512     lastfilestep = 0;
 513     laststep     = 0;
 514
 515     auto files = gmx::copyOf(opt2fns("-f", NFILE, fnm));
 516
 517     if (files.empty())
 518     {
 519         gmx_fatal(FARGS, "No input files!");
 520     }
 521
 522     snew(settime, files.size() + 1);
 523     snew(readtime, files.size() + 1);
 524     snew(cont_type, files.size() + 1);
 525
 526     nre = scan_ene_files(files, readtime, &timestep, &nremax);
 527     edit_files(files, readtime, settime, cont_type, bSetTime, bSort);
 528
 529     ee_sum_nsteps = 0;
 530     ee_sum_nsum   = 0;
 531     snew(ee_sum, nremax);
 532
 533     snew(fr, 1);
 534     snew(fro, 1);
 535     fro->t   = -1e20;
 536     fro->nre = nre;
 537     snew(fro->ener, nremax);
 538
 539     noutfr = 0;
 540     bFirst = TRUE;
 541
 542     last_t = fro->t;
 543     for (size_t f = 0; f < files.size(); f++)
 544     {
 545         bNewFile   = TRUE;
 546         bNewOutput = TRUE;
 547         in         = open_enx(files[f].c_str(), "r");
 548         enm        = nullptr;
 549         do_enxnms(in, &this_nre, &enm);
 550         if (f == 0)
 551         {
 552             if (scalefac != 1)
 553             {
 554                 set = select_it(nre, enm, &nset);
 555             }
 556
 557             /* write names to the output file */
 558             out = open_enx(opt2fn("-o", NFILE, fnm), "w");
 559             do_enxnms(out, &nre, &enm);
 560         }
 561
 562         /* start reading from the next file */
 563         while ((fro->t <= (settime[f+1] + GMX_REAL_EPS)) &&
 564                do_enx(in, fr))
 565         {
 566             if (bNewFile)
 567             {
 568                 startstep_file = fr->step;
 569                 tadjust        = settime[f] - fr->t;
 570                 if (cont_type[f+1] == TIME_LAST)
 571                 {
 572                     settime[f+1]   = readtime[f+1]-readtime[f]+settime[f];
 573                     cont_type[f+1] = TIME_EXPLICIT;
 574                 }
 575                 bNewFile = FALSE;
 576             }
 577
 578             if (tadjust + fr->t <= last_t)
 579             {
 580                 /* Skip this frame, since we already have it / past it */
 581                 if (debug)
 582                 {
 583                     fprintf(debug, "fr->step %s, fr->t %.4f\n",
 584                             gmx_step_str(fr->step, buf), fr->t);
 585                     fprintf(debug, "tadjust %12.6e + fr->t %12.6e <= t %12.6e\n",
 586                             tadjust, fr->t, last_t);
 587                 }
 588                 continue;
 589             }
 590
 591             fro->step = lastfilestep + fr->step - startstep_file;
 592             fro->t    = tadjust  + fr->t;
 593
 594             bWrite = ((begin < 0 || (fro->t >= begin-GMX_REAL_EPS)) &&
 595                       (end  < 0  || (fro->t <= end  +GMX_REAL_EPS)) &&
 596                       (fro->t <= settime[f+1]+0.5*timestep));
 597
 598             if (debug)
 599             {
 600                 fprintf(debug,
 601                         "fr->step %s, fr->t %.4f, fro->step %s fro->t %.4f, w %s\n",
 602                         gmx_step_str(fr->step, buf), fr->t,
 603                         gmx_step_str(fro->step, buf2), fro->t, gmx::boolToString(bWrite));
 604             }
 605
 606             if (bError)
 607             {
 608                 if ((end > 0) && (fro->t > end+GMX_REAL_EPS))
 609                 {
 610                     f = files.size();
 611                     break;
 612                 }
 613             }
 614
 615             if (fro->t >= begin-GMX_REAL_EPS)
 616             {
 617                 if (bFirst)
 618                 {
 619                     bFirst    = FALSE;
 620                 }
 621                 if (bWrite)
 622                 {
 623                     update_ee_sum(nre, &ee_sum_step, &ee_sum_nsteps, &ee_sum_nsum, ee_sum,
 624                                   fr, fro->step);
 625                 }
 626             }
 627
 628             /* determine if we should write it */
 629             if (bWrite && (delta_t == 0 || bRmod(fro->t, toffset, delta_t)))
 630             {
 631                 laststep = fro->step;
 632                 last_t   = fro->t;
 633                 if (bNewOutput)
 634                 {
 635                     bNewOutput = FALSE;
 636                     fprintf(stderr, "\nContinue writing frames from t=%g, step=%s\n",
 637                             fro->t, gmx_step_str(fro->step, buf));
 638                 }
 639
 640                 /* Copy the energies */
 641                 for (i = 0; i < nre; i++)
 642                 {
 643                     fro->ener[i].e = fr->ener[i].e;
 644                 }
 645
 646                 fro->nsteps = ee_sum_nsteps;
 647                 fro->dt     = fr->dt;
 648
 649                 if (ee_sum_nsum <= 1)
 650                 {
 651                     fro->nsum = 0;
 652                 }
 653                 else
 654                 {
 655                     fro->nsum = int64_to_int(ee_sum_nsum,
 656                                              "energy average summation");
 657                     /* Copy the energy sums */
 658                     for (i = 0; i < nre; i++)
 659                     {
 660                         fro->ener[i].esum = ee_sum[i].esum;
 661                         fro->ener[i].eav  = ee_sum[i].eav;
 662                     }
 663                 }
 664                 /* We wrote the energies, so reset the counts */
 665                 ee_sum_nsteps = 0;
 666                 ee_sum_nsum   = 0;
 667
 668                 if (scalefac != 1)
 669                 {
 670                     for (kkk = 0; kkk < nset; kkk++)
 671                     {
 672                         fro->ener[set[kkk]].e    *= scalefac;
 673                         if (fro->nsum > 0)
 674                         {
 675                             fro->ener[set[kkk]].eav  *= scalefac*scalefac;
 676                             fro->ener[set[kkk]].esum *= scalefac;
 677                         }
 678                     }
 679                 }
 680                 /* Copy restraint stuff */
 681                 /*fro->ndisre       = fr->ndisre;
 682                    fro->disre_rm3tav = fr->disre_rm3tav;
 683                    fro->disre_rt     = fr->disre_rt;*/
 684                 fro->nblock       = fr->nblock;
 685                 /*fro->nr           = fr->nr;*/
 686                 fro->block        = fr->block;
 687
 688                 /* check if we have blocks with delta_h data and are throwing
 689                    away data */
 690                 if (fro->nblock > 0)
 691                 {
 692                     if (remove_dh)
 693                     {
 694                         int i;
 695                         if (!blocks || nblocks_alloc < fr->nblock)
 696                         {
 697                             /* we pre-allocate the blocks */
 698                             nblocks_alloc = fr->nblock;
 699                             snew(blocks, nblocks_alloc);
 700                         }
 701                         nblocks = 0; /* number of blocks so far */
 702
 703                         for (i = 0; i < fr->nblock; i++)
 704                         {
 705                             if ( (fr->block[i].id != enxDHCOLL) &&
 706                                  (fr->block[i].id != enxDH) &&
 707                                  (fr->block[i].id != enxDHHIST) )
 708                             {
 709                                 /* copy everything verbatim */
 710                                 blocks[nblocks] = fr->block[i];
 711                                 nblocks++;
 712                             }
 713                         }
 714                         /* now set the block pointer to the new blocks */
 715                         fro->nblock = nblocks;
 716                         fro->block  = blocks;
 717                     }
 718                     else if (delta_t > 0)
 719                     {
 720                         if (!warned_about_dh)
 721                         {
 722                             for (i = 0; i < fr->nblock; i++)
 723                             {
 724                                 if (fr->block[i].id == enxDH ||
 725                                     fr->block[i].id == enxDHHIST)
 726                                 {
 727                                     int size;
 728                                     if (fr->block[i].id == enxDH)
 729                                     {
 730                                         size = fr->block[i].sub[2].nr;
 731                                     }
 732                                     else
 733                                     {
 734                                         size = fr->nsteps;
 735                                     }
 736                                     if (size > 0)
 737                                     {
 738                                         printf("\nWARNING: %s contains delta H blocks or histograms for which\n"
 739                                                "         some data is thrown away on a block-by-block basis, where each block\n"
 740                                                "         contains up to %d samples.\n"
 741                                                "         This is almost certainly not what you want.\n"
 742                                                "         Use the -rmdh option to throw all delta H samples away.\n"
 743                                                "         Use g_energy -odh option to extract these samples.\n",
 744                                                files[f].c_str(), size);
 745                                         warned_about_dh = TRUE;
 746                                         break;
 747                                     }
 748                                 }
 749                             }
 750                         }
 751                     }
 752                 }
 753
 754                 do_enx(out, fro);
 755                 if (noutfr % 1000 == 0)
 756                 {
 757                     fprintf(stderr, "Writing frame time %g    ", fro->t);
 758                 }
 759                 noutfr++;
 760             }
 761         }
 762         if (f == files.size())
 763         {
 764             f--;
 765         }
 766         printf("\nLast step written from %s: t %g, step %s\n",
 767                files[f].c_str(), last_t, gmx_step_str(laststep, buf));
 768         lastfilestep = laststep;
 769
 770         /* set the next time from the last in previous file */
 771         if (cont_type[f+1] == TIME_CONTINUE)
 772         {
 773             settime[f+1] = fro->t;
 774             /* in this case we have already written the last frame of
 775              * previous file, so update begin to avoid doubling it
 776              * with the start of the next file
 777              */
 778             begin = fro->t+0.5*timestep;
 779             /* cont_type[f+1]==TIME_EXPLICIT; */
 780         }
 781
 782         if ((fro->t < end) && (f < files.size() - 1) &&
 783             (fro->t < settime[f+1]-1.5*timestep))
 784         {
 785             fprintf(stderr,
 786                     "\nWARNING: There might be a gap around t=%g\n", fro->t);
 787         }
 788
 789         /* move energies to lastee */
 790         close_enx(in);
 791         free_enxnms(this_nre, enm);
 792
 793         fprintf(stderr, "\n");
 794     }
 795     if (noutfr == 0)
 796     {
 797         fprintf(stderr, "No frames written.\n");
 798     }
 799     else
 800     {
 801         fprintf(stderr, "Last frame written was at step %s, time %f\n",
 802                 gmx_step_str(fro->step, buf), fro->t);
 803         fprintf(stderr, "Wrote %d frames\n", noutfr);
 804     }
 805
 806     return 0;
 807 }