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,2021, 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(),
 191                         nresav,
 192                         files[f].c_str(),
 193                         fr->nre);
 194                 fprintf(stderr,
 195                         "\nContinue conversion using only the first %d terms (n/y)?\n"
 196                         "(you should be sure that the energy terms match)\n",
 197                         nremin);
 198                 if (nullptr == fgets(inputstring, STRLEN - 1, stdin))
 199                 {
 200                     gmx_fatal(FARGS, "Error reading user input");
 201                 }
 202                 if (inputstring[0] != 'y' && inputstring[0] != 'Y')
 203                 {
 204                     fprintf(stderr, "Will not convert\n");
 205                     exit(0);
 206                 }
 207                 nresav = fr->nre;
 208             }
 209             do_enx(in, fr);
 210             readtime[f] = fr->t;
 211             close_enx(in);
 212         }
 213         fprintf(stderr, "\n");
 214         free_enxnms(nre, enm);
 215     }
 216
 217     free_enxframe(fr);
 218     sfree(fr);
 219
 220     return nremin;
 221 }
 222
 223
 224 static void edit_files(gmx::ArrayRef<std::string> files,
 225                        real*                      readtime,
 226                        real*                      settime,
 227                        int*                       cont_type,
 228                        gmx_bool                   bSetTime,
 229                        gmx_bool                   bSort)
 230 {
 231     gmx_bool ok;
 232     char     inputstring[STRLEN], *chptr;
 233
 234     if (bSetTime)
 235     {
 236         if (files.size() == 1)
 237         {
 238             fprintf(stderr, "\n\nEnter the new start time:\n\n");
 239         }
 240         else
 241         {
 242             fprintf(stderr,
 243                     "\n\nEnter the new start time for each file.\n"
 244                     "There are two special options, both disables sorting:\n\n"
 245                     "c (continue) - The start time is taken from the end\n"
 246                     "of the previous file. Use it when your continuation run\n"
 247                     "restarts with t=0 and there is no overlap.\n\n"
 248                     "l (last) - The time in this file will be changed the\n"
 249                     "same amount as in the previous. Use it when the time in the\n"
 250                     "new run continues from the end of the previous one,\n"
 251                     "since this takes possible overlap into account.\n\n");
 252         }
 253
 254         fprintf(stderr,
 255                 "          File             Current start       New start\n"
 256                 "---------------------------------------------------------\n");
 257
 258         for (gmx::index i = 0; i < files.ssize(); i++)
 259         {
 260             fprintf(stderr, "%25s   %10.3f             ", files[i].c_str(), readtime[i]);
 261             ok = FALSE;
 262             do
 263             {
 264                 if (nullptr == fgets(inputstring, STRLEN - 1, stdin))
 265                 {
 266                     gmx_fatal(FARGS, "Error reading user input");
 267                 }
 268                 inputstring[std::strlen(inputstring) - 1] = 0;
 269
 270                 if (inputstring[0] == 'c' || inputstring[0] == 'C')
 271                 {
 272                     cont_type[i] = TIME_CONTINUE;
 273                     bSort        = FALSE;
 274                     ok           = TRUE;
 275                     settime[i]   = FLT_MAX;
 276                 }
 277                 else if (inputstring[0] == 'l' || inputstring[0] == 'L')
 278                 {
 279                     cont_type[i] = TIME_LAST;
 280                     bSort        = FALSE;
 281                     ok           = TRUE;
 282                     settime[i]   = FLT_MAX;
 283                 }
 284                 else
 285                 {
 286                     settime[i] = strtod(inputstring, &chptr);
 287                     if (chptr == inputstring)
 288                     {
 289                         fprintf(stderr, "Try that again: ");
 290                     }
 291                     else
 292                     {
 293                         cont_type[i] = TIME_EXPLICIT;
 294                         ok           = TRUE;
 295                     }
 296                 }
 297             } while (!ok);
 298         }
 299         if (cont_type[0] != TIME_EXPLICIT)
 300         {
 301             cont_type[0] = TIME_EXPLICIT;
 302             settime[0]   = 0;
 303         }
 304     }
 305     else
 306     {
 307         for (gmx::index i = 0; i < files.ssize(); i++)
 308         {
 309             settime[i] = readtime[i];
 310         }
 311     }
 312
 313     if (bSort && files.size() > 1)
 314     {
 315         sort_files(files, settime);
 316     }
 317     else
 318     {
 319         fprintf(stderr, "Sorting disabled.\n");
 320     }
 321
 322
 323     /* Write out the new order and start times */
 324     fprintf(stderr,
 325             "\nSummary of files and start times used:\n\n"
 326             "          File                Start time\n"
 327             "-----------------------------------------\n");
 328     for (gmx::index i = 0; i < files.ssize(); i++)
 329     {
 330         switch (cont_type[i])
 331         {
 332             case TIME_EXPLICIT:
 333                 fprintf(stderr, "%25s   %10.3f\n", files[i].c_str(), settime[i]);
 334                 break;
 335             case TIME_CONTINUE:
 336                 fprintf(stderr, "%25s        Continue from end of last file\n", files[i].c_str());
 337                 break;
 338             case TIME_LAST:
 339                 fprintf(stderr, "%25s        Change by same amount as last file\n", files[i].c_str());
 340                 break;
 341         }
 342     }
 343     fprintf(stderr, "\n");
 344
 345     settime[files.size()]   = FLT_MAX;
 346     cont_type[files.size()] = TIME_EXPLICIT;
 347     readtime[files.size()]  = FLT_MAX;
 348 }
 349
 350
 351 static void update_ee_sum(int         nre,
 352                           int64_t*    ee_sum_step,
 353                           int64_t*    ee_sum_nsteps,
 354                           int64_t*    ee_sum_nsum,
 355                           t_energy*   ee_sum,
 356                           t_enxframe* fr,
 357                           int         out_step)
 358 {
 359     int64_t nsteps, nsum, fr_nsum;
 360     int     i;
 361
 362     nsteps = *ee_sum_nsteps;
 363     nsum   = *ee_sum_nsum;
 364
 365     fr_nsum = fr->nsum;
 366     if (fr_nsum == 0)
 367     {
 368         fr_nsum = 1;
 369     }
 370
 371     if (nsteps == 0)
 372     {
 373         if (fr_nsum == 1)
 374         {
 375             for (i = 0; i < nre; i++)
 376             {
 377                 ee_sum[i].esum = fr->ener[i].e;
 378                 ee_sum[i].eav  = 0;
 379             }
 380         }
 381         else
 382         {
 383             for (i = 0; i < nre; i++)
 384             {
 385                 ee_sum[i].esum = fr->ener[i].esum;
 386                 ee_sum[i].eav  = fr->ener[i].eav;
 387             }
 388         }
 389         nsteps = fr->nsteps;
 390         nsum   = fr_nsum;
 391     }
 392     else if (out_step + *ee_sum_nsum - *ee_sum_step == nsteps + fr->nsteps)
 393     {
 394         if (fr_nsum == 1)
 395         {
 396             for (i = 0; i < nre; i++)
 397             {
 398                 ee_sum[i].eav += gmx::square(ee_sum[i].esum / nsum
 399                                              - (ee_sum[i].esum + fr->ener[i].e) / (nsum + 1))
 400                                  * nsum * (nsum + 1);
 401                 ee_sum[i].esum += fr->ener[i].e;
 402             }
 403         }
 404         else
 405         {
 406             for (i = 0; i < fr->nre; i++)
 407             {
 408                 ee_sum[i].eav += fr->ener[i].eav
 409                                  + gmx::square(ee_sum[i].esum / nsum
 410                                                - (ee_sum[i].esum + fr->ener[i].esum) / (nsum + fr->nsum))
 411                                            * nsum * (nsum + fr->nsum) / static_cast<double>(fr->nsum);
 412                 ee_sum[i].esum += fr->ener[i].esum;
 413             }
 414         }
 415         nsteps += fr->nsteps;
 416         nsum += fr_nsum;
 417     }
 418     else
 419     {
 420         if (fr->nsum != 0)
 421         {
 422             fprintf(stderr, "\nWARNING: missing energy sums at time %f\n", fr->t);
 423         }
 424         nsteps = 0;
 425         nsum   = 0;
 426     }
 427
 428     *ee_sum_step   = out_step;
 429     *ee_sum_nsteps = nsteps;
 430     *ee_sum_nsum   = nsum;
 431 }
 432
 433 int gmx_eneconv(int argc, char* argv[])
 434 {
 435     const char* desc[] = {
 436         "With [IT]multiple files[it] specified for the [TT]-f[tt] option:[PAR]",
 437         "Concatenates several energy files in sorted order.",
 438         "In the case of double time frames, the one",
 439         "in the later file is used. By specifying [TT]-settime[tt] you will be",
 440         "asked for the start time of each file. The input files are taken",
 441         "from the command line,",
 442         "such that the command [TT]gmx eneconv -f *.edr -o fixed.edr[tt] should do",
 443         "the trick. [PAR]",
 444         "With [IT]one file[it] specified for [TT]-f[tt]:[PAR]",
 445         "Reads one energy file and writes another, applying the [TT]-dt[tt],",
 446         "[TT]-offset[tt], [TT]-t0[tt] and [TT]-settime[tt] options and",
 447         "converting to a different format if necessary (indicated by file",
 448         "extensions).[PAR]",
 449         "[TT]-settime[tt] is applied first, then [TT]-dt[tt]/[TT]-offset[tt]",
 450         "followed by [TT]-b[tt] and [TT]-e[tt] to select which frames to write."
 451     };
 452     const char* bugs[] = {
 453         "When combining trajectories the sigma and E^2 (necessary for statistics) are not "
 454         "updated correctly. Only the actual energy is correct. One thus has to compute "
 455         "statistics in another way."
 456     };
 457     ener_file_t  in = nullptr, out = nullptr;
 458     gmx_enxnm_t* enm = nullptr;
 459 #if 0
 460     ener_file_t       in, out = NULL;
 461     gmx_enxnm_t      *enm = NULL;
 462 #endif
 463     t_enxframe *      fr, *fro;
 464     int64_t           ee_sum_step = 0, ee_sum_nsteps, ee_sum_nsum;
 465     t_energy*         ee_sum;
 466     int64_t           lastfilestep, laststep, startstep_file = 0;
 467     int               noutfr;
 468     int               nre, nremax, this_nre, i, kkk, nset, *set = nullptr;
 469     double            last_t;
 470     real *            readtime, *settime, timestep, tadjust;
 471     char              buf[22], buf2[22];
 472     int*              cont_type;
 473     gmx_bool          bNewFile, bFirst, bNewOutput;
 474     gmx_output_env_t* oenv;
 475     gmx_bool          warned_about_dh = FALSE;
 476     t_enxblock*       blocks          = nullptr;
 477     int               nblocks         = 0;
 478     int               nblocks_alloc   = 0;
 479
 480     t_filenm fnm[] = {
 481         { efEDR, "-f", nullptr, ffRDMULT },
 482         { efEDR, "-o", "fixed", ffWRITE },
 483     };
 484
 485 #define NFILE asize(fnm)
 486     gmx_bool        bWrite;
 487     static real     delta_t = 0.0, toffset = 0, scalefac = 1;
 488     static gmx_bool bSetTime = FALSE;
 489     static gmx_bool bSort = TRUE, bError = TRUE;
 490     static real     begin     = -1;
 491     static real     end       = -1;
 492     gmx_bool        remove_dh = FALSE;
 493
 494     t_pargs pa[] = {
 495         { "-b", FALSE, etREAL, { &begin }, "First time to use" },
 496         { "-e", FALSE, etREAL, { &end }, "Last time to use" },
 497         { "-dt", FALSE, etREAL, { &delta_t }, "Only write out frame when t MOD dt = offset" },
 498         { "-offset", FALSE, etREAL, { &toffset }, "Time offset for [TT]-dt[tt] option" },
 499         { "-settime", FALSE, etBOOL, { &bSetTime }, "Change starting time interactively" },
 500         { "-sort", FALSE, etBOOL, { &bSort }, "Sort energy files (not frames)" },
 501         { "-rmdh", FALSE, etBOOL, { &remove_dh }, "Remove free energy block data" },
 502         { "-scalefac", FALSE, etREAL, { &scalefac }, "Multiply energy component by this factor" },
 503         { "-error", FALSE, etBOOL, { &bError }, "Stop on errors in the file" }
 504     };
 505
 506     if (!parse_common_args(
 507                 &argc, argv, 0, NFILE, fnm, asize(pa), pa, asize(desc), desc, asize(bugs), bugs, &oenv))
 508     {
 509         return 0;
 510     }
 511     fprintf(stdout,
 512             "Note that major changes are planned in future for "
 513             "eneconv, to improve usability and utility.");
 514
 515     tadjust      = 0;
 516     nremax       = 0;
 517     nset         = 0;
 518     timestep     = 0.0;
 519     lastfilestep = 0;
 520     laststep     = 0;
 521
 522     auto files = gmx::copyOf(opt2fns("-f", NFILE, fnm));
 523
 524     if (files.empty())
 525     {
 526         gmx_fatal(FARGS, "No input files!");
 527     }
 528
 529     snew(settime, files.size() + 1);
 530     snew(readtime, files.size() + 1);
 531     snew(cont_type, files.size() + 1);
 532
 533     nre = scan_ene_files(files, readtime, &timestep, &nremax);
 534     edit_files(files, readtime, settime, cont_type, bSetTime, bSort);
 535
 536     ee_sum_nsteps = 0;
 537     ee_sum_nsum   = 0;
 538     snew(ee_sum, nremax);
 539
 540     snew(fr, 1);
 541     snew(fro, 1);
 542     fro->t   = -1e20;
 543     fro->nre = nre;
 544     snew(fro->ener, nremax);
 545
 546     noutfr = 0;
 547     bFirst = TRUE;
 548
 549     last_t = fro->t;
 550     for (size_t f = 0; f < files.size(); f++)
 551     {
 552         bNewFile   = TRUE;
 553         bNewOutput = TRUE;
 554         in         = open_enx(files[f].c_str(), "r");
 555         enm        = nullptr;
 556         do_enxnms(in, &this_nre, &enm);
 557         if (f == 0)
 558         {
 559             if (scalefac != 1)
 560             {
 561                 set = select_it(nre, enm, &nset);
 562             }
 563
 564             /* write names to the output file */
 565             out = open_enx(opt2fn("-o", NFILE, fnm), "w");
 566             do_enxnms(out, &nre, &enm);
 567         }
 568
 569         /* start reading from the next file */
 570         while ((fro->t <= (settime[f + 1] + GMX_REAL_EPS)) && do_enx(in, fr))
 571         {
 572             if (bNewFile)
 573             {
 574                 startstep_file = fr->step;
 575                 tadjust        = settime[f] - fr->t;
 576                 if (cont_type[f + 1] == TIME_LAST)
 577                 {
 578                     settime[f + 1]   = readtime[f + 1] - readtime[f] + settime[f];
 579                     cont_type[f + 1] = TIME_EXPLICIT;
 580                 }
 581                 bNewFile = FALSE;
 582             }
 583
 584             if (tadjust + fr->t <= last_t)
 585             {
 586                 /* Skip this frame, since we already have it / past it */
 587                 if (debug)
 588                 {
 589                     fprintf(debug, "fr->step %s, fr->t %.4f\n", gmx_step_str(fr->step, buf), fr->t);
 590                     fprintf(debug, "tadjust %12.6e + fr->t %12.6e <= t %12.6e\n", tadjust, fr->t, last_t);
 591                 }
 592                 continue;
 593             }
 594
 595             fro->step = lastfilestep + fr->step - startstep_file;
 596             fro->t    = tadjust + fr->t;
 597
 598             bWrite = ((begin < 0 || (fro->t >= begin - GMX_REAL_EPS))
 599                       && (end < 0 || (fro->t <= end + GMX_REAL_EPS))
 600                       && (fro->t <= settime[f + 1] + 0.5 * timestep));
 601
 602             if (debug)
 603             {
 604                 fprintf(debug,
 605                         "fr->step %s, fr->t %.4f, fro->step %s fro->t %.4f, w %s\n",
 606                         gmx_step_str(fr->step, buf),
 607                         fr->t,
 608                         gmx_step_str(fro->step, buf2),
 609                         fro->t,
 610                         gmx::boolToString(bWrite));
 611             }
 612
 613             if (bError)
 614             {
 615                 if ((end > 0) && (fro->t > end + GMX_REAL_EPS))
 616                 {
 617                     f = files.size();
 618                     break;
 619                 }
 620             }
 621
 622             if (fro->t >= begin - GMX_REAL_EPS)
 623             {
 624                 if (bFirst)
 625                 {
 626                     bFirst = FALSE;
 627                 }
 628                 if (bWrite)
 629                 {
 630                     update_ee_sum(nre, &ee_sum_step, &ee_sum_nsteps, &ee_sum_nsum, ee_sum, fr, fro->step);
 631                 }
 632             }
 633
 634             /* determine if we should write it */
 635             if (bWrite && (delta_t == 0 || bRmod(fro->t, toffset, delta_t)))
 636             {
 637                 laststep = fro->step;
 638                 last_t   = fro->t;
 639                 if (bNewOutput)
 640                 {
 641                     bNewOutput = FALSE;
 642                     fprintf(stderr,
 643                             "\nContinue writing frames from t=%g, step=%s\n",
 644                             fro->t,
 645                             gmx_step_str(fro->step, buf));
 646                 }
 647
 648                 /* Copy the energies */
 649                 for (i = 0; i < nre; i++)
 650                 {
 651                     fro->ener[i].e = fr->ener[i].e;
 652                 }
 653
 654                 fro->nsteps = ee_sum_nsteps;
 655                 fro->dt     = fr->dt;
 656
 657                 if (ee_sum_nsum <= 1)
 658                 {
 659                     fro->nsum = 0;
 660                 }
 661                 else
 662                 {
 663                     fro->nsum = int64_to_int(ee_sum_nsum, "energy average summation");
 664                     /* Copy the energy sums */
 665                     for (i = 0; i < nre; i++)
 666                     {
 667                         fro->ener[i].esum = ee_sum[i].esum;
 668                         fro->ener[i].eav  = ee_sum[i].eav;
 669                     }
 670                 }
 671                 /* We wrote the energies, so reset the counts */
 672                 ee_sum_nsteps = 0;
 673                 ee_sum_nsum   = 0;
 674
 675                 if (scalefac != 1)
 676                 {
 677                     for (kkk = 0; kkk < nset; kkk++)
 678                     {
 679                         fro->ener[set[kkk]].e *= scalefac;
 680                         if (fro->nsum > 0)
 681                         {
 682                             fro->ener[set[kkk]].eav *= scalefac * scalefac;
 683                             fro->ener[set[kkk]].esum *= scalefac;
 684                         }
 685                     }
 686                 }
 687                 /* Copy restraint stuff */
 688                 /*fro->ndisre       = fr->ndisre;
 689                    fro->disre_rm3tav = fr->disre_rm3tav;
 690                    fro->disre_rt     = fr->disre_rt;*/
 691                 fro->nblock = fr->nblock;
 692                 /*fro->nr           = fr->nr;*/
 693                 fro->block = fr->block;
 694
 695                 /* check if we have blocks with delta_h data and are throwing
 696                    away data */
 697                 if (fro->nblock > 0)
 698                 {
 699                     if (remove_dh)
 700                     {
 701                         int i;
 702                         if (!blocks || nblocks_alloc < fr->nblock)
 703                         {
 704                             /* we pre-allocate the blocks */
 705                             nblocks_alloc = fr->nblock;
 706                             snew(blocks, nblocks_alloc);
 707                         }
 708                         nblocks = 0; /* number of blocks so far */
 709
 710                         for (i = 0; i < fr->nblock; i++)
 711                         {
 712                             if ((fr->block[i].id != enxDHCOLL) && (fr->block[i].id != enxDH)
 713                                 && (fr->block[i].id != enxDHHIST))
 714                             {
 715                                 /* copy everything verbatim */
 716                                 blocks[nblocks] = fr->block[i];
 717                                 nblocks++;
 718                             }
 719                         }
 720                         /* now set the block pointer to the new blocks */
 721                         fro->nblock = nblocks;
 722                         fro->block  = blocks;
 723                     }
 724                     else if (delta_t > 0)
 725                     {
 726                         if (!warned_about_dh)
 727                         {
 728                             for (i = 0; i < fr->nblock; i++)
 729                             {
 730                                 if (fr->block[i].id == enxDH || fr->block[i].id == enxDHHIST)
 731                                 {
 732                                     int size;
 733                                     if (fr->block[i].id == enxDH)
 734                                     {
 735                                         size = fr->block[i].sub[2].nr;
 736                                     }
 737                                     else
 738                                     {
 739                                         size = fr->nsteps;
 740                                     }
 741                                     if (size > 0)
 742                                     {
 743                                         printf("\nWARNING: %s contains delta H blocks or "
 744                                                "histograms for which\n"
 745                                                "         some data is thrown away on a "
 746                                                "block-by-block basis, where each block\n"
 747                                                "         contains up to %d samples.\n"
 748                                                "         This is almost certainly not what you "
 749                                                "want.\n"
 750                                                "         Use the -rmdh option to throw all delta H "
 751                                                "samples away.\n"
 752                                                "         Use g_energy -odh option to extract these "
 753                                                "samples.\n",
 754                                                files[f].c_str(),
 755                                                size);
 756                                         warned_about_dh = TRUE;
 757                                         break;
 758                                     }
 759                                 }
 760                             }
 761                         }
 762                     }
 763                 }
 764
 765                 do_enx(out, fro);
 766                 if (noutfr % 1000 == 0)
 767                 {
 768                     fprintf(stderr, "Writing frame time %g    ", fro->t);
 769                 }
 770                 noutfr++;
 771             }
 772         }
 773         if (f == files.size())
 774         {
 775             f--;
 776         }
 777         printf("\nLast step written from %s: t %g, step %s\n",
 778                files[f].c_str(),
 779                last_t,
 780                gmx_step_str(laststep, buf));
 781         lastfilestep = laststep;
 782
 783         /* set the next time from the last in previous file */
 784         if (cont_type[f + 1] == TIME_CONTINUE)
 785         {
 786             settime[f + 1] = fro->t;
 787             /* in this case we have already written the last frame of
 788              * previous file, so update begin to avoid doubling it
 789              * with the start of the next file
 790              */
 791             begin = fro->t + 0.5 * timestep;
 792             /* cont_type[f+1]==TIME_EXPLICIT; */
 793         }
 794
 795         if ((fro->t < end) && (f < files.size() - 1) && (fro->t < settime[f + 1] - 1.5 * timestep))
 796         {
 797             fprintf(stderr, "\nWARNING: There might be a gap around t=%g\n", fro->t);
 798         }
 799
 800         /* move energies to lastee */
 801         close_enx(in);
 802         free_enxnms(this_nre, enm);
 803
 804         fprintf(stderr, "\n");
 805     }
 806     if (noutfr == 0)
 807     {
 808         fprintf(stderr, "No frames written.\n");
 809     }
 810     else
 811     {
 812         fprintf(stderr, "Last frame written was at step %s, time %f\n", gmx_step_str(fro->step, buf), fro->t);
 813         fprintf(stderr, "Wrote %d frames\n", noutfr);
 814     }
 815
 816     return 0;
 817 }