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