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