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