src/gromacs/gmxana/gmx_eneconv.cpp

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