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, 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/utility/arraysize.h"
  54 #include "gromacs/utility/cstringutil.h"
  55 #include "gromacs/utility/fatalerror.h"
  56 #include "gromacs/utility/int64_to_int.h"
  57 #include "gromacs/utility/smalloc.h"
  58
  59 #define TIME_EXPLICIT 0
  60 #define TIME_CONTINUE 1
  61 #define TIME_LAST     2
  62 #ifndef FLT_MAX
  63 #define FLT_MAX 1e36
  64 #endif
  65
  66 static int *select_it(int nre, gmx_enxnm_t *nm, int *nset)
  67 {
  68     gmx_bool *bE;
  69     int       n, k, j, i;
  70     int      *set;
  71     gmx_bool  bVerbose = TRUE;
  72
  73     if ((getenv("GMX_ENER_VERBOSE")) != NULL)
  74     {
  75         bVerbose = FALSE;
  76     }
  77
  78     fprintf(stderr, "Select the terms you want to scale from the following list\n");
  79     fprintf(stderr, "End your selection with 0\n");
  80
  81     if (bVerbose)
  82     {
  83         for (k = 0; (k < nre); )
  84         {
  85             for (j = 0; (j < 4) && (k < nre); j++, k++)
  86             {
  87                 fprintf(stderr, " %3d=%14s", k+1, nm[k].name);
  88             }
  89             fprintf(stderr, "\n");
  90         }
  91     }
  92
  93     snew(bE, nre);
  94     do
  95     {
  96         if (1 != scanf("%d", &n))
  97         {
  98             gmx_fatal(FARGS, "Cannot read energy term");
  99         }
 100         if ((n > 0) && (n <= nre))
 101         {
 102             bE[n-1] = TRUE;
 103         }
 104     }
 105     while (n != 0);
 106
 107     snew(set, nre);
 108     for (i = (*nset) = 0; (i < nre); i++)
 109     {
 110         if (bE[i])
 111         {
 112             set[(*nset)++] = i;
 113         }
 114     }
 115
 116     sfree(bE);
 117
 118     return set;
 119 }
 120
 121 static void sort_files(char **fnms, real *settime, int nfile)
 122 {
 123     int   i, j, minidx;
 124     real  timeswap;
 125     char *chptr;
 126
 127     for (i = 0; i < nfile; i++)
 128     {
 129         minidx = i;
 130         for (j = i+1; j < nfile; j++)
 131         {
 132             if (settime[j] < settime[minidx])
 133             {
 134                 minidx = j;
 135             }
 136         }
 137         if (minidx != i)
 138         {
 139             timeswap        = settime[i];
 140             settime[i]      = settime[minidx];
 141             settime[minidx] = timeswap;
 142             chptr           = fnms[i];
 143             fnms[i]         = fnms[minidx];
 144             fnms[minidx]    = chptr;
 145         }
 146     }
 147 }
 148
 149
 150 static int scan_ene_files(char **fnms, int nfiles,
 151                           real *readtime, real *timestep, int *nremax)
 152 {
 153     /* Check number of energy terms and start time of all files */
 154     int          f, 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 (f = 0; f < nfiles; f++)
 164     {
 165         in  = open_enx(fnms[f], "r");
 166         enm = NULL;
 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", fnms[f-1], nresav, fnms[f], 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 (NULL == 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(char **fnms, int nfiles, real *readtime,
 221                        real *settime, int *cont_type, gmx_bool bSetTime, gmx_bool bSort)
 222 {
 223     int      i;
 224     gmx_bool ok;
 225     char     inputstring[STRLEN], *chptr;
 226
 227     if (bSetTime)
 228     {
 229         if (nfiles == 1)
 230         {
 231             fprintf(stderr, "\n\nEnter the new start time:\n\n");
 232         }
 233         else
 234         {
 235             fprintf(stderr, "\n\nEnter the new start time for each file.\n"
 236                     "There are two special options, both disables sorting:\n\n"
 237                     "c (continue) - The start time is taken from the end\n"
 238                     "of the previous file. Use it when your continuation run\n"
 239                     "restarts with t=0 and there is no overlap.\n\n"
 240                     "l (last) - The time in this file will be changed the\n"
 241                     "same amount as in the previous. Use it when the time in the\n"
 242                     "new run continues from the end of the previous one,\n"
 243                     "since this takes possible overlap into account.\n\n");
 244         }
 245
 246         fprintf(stderr, "          File             Current start       New start\n"
 247                 "---------------------------------------------------------\n");
 248
 249         for (i = 0; i < nfiles; i++)
 250         {
 251             fprintf(stderr, "%25s   %10.3f             ", fnms[i], readtime[i]);
 252             ok = FALSE;
 253             do
 254             {
 255                 if (NULL == fgets(inputstring, STRLEN-1, stdin))
 256                 {
 257                     gmx_fatal(FARGS, "Error reading user input");
 258                 }
 259                 inputstring[std::strlen(inputstring)-1] = 0;
 260
 261                 if (inputstring[0] == 'c' || inputstring[0] == 'C')
 262                 {
 263                     cont_type[i] = TIME_CONTINUE;
 264                     bSort        = FALSE;
 265                     ok           = TRUE;
 266                     settime[i]   = FLT_MAX;
 267                 }
 268                 else if (inputstring[0] == 'l' ||
 269                          inputstring[0] == 'L')
 270                 {
 271                     cont_type[i] = TIME_LAST;
 272                     bSort        = FALSE;
 273                     ok           = TRUE;
 274                     settime[i]   = FLT_MAX;
 275                 }
 276                 else
 277                 {
 278                     settime[i] = strtod(inputstring, &chptr);
 279                     if (chptr == inputstring)
 280                     {
 281                         fprintf(stderr, "Try that again: ");
 282                     }
 283                     else
 284                     {
 285                         cont_type[i] = TIME_EXPLICIT;
 286                         ok           = TRUE;
 287                     }
 288                 }
 289             }
 290             while (!ok);
 291         }
 292         if (cont_type[0] != TIME_EXPLICIT)
 293         {
 294             cont_type[0] = TIME_EXPLICIT;
 295             settime[0]   = 0;
 296         }
 297     }
 298     else
 299     {
 300         for (i = 0; i < nfiles; i++)
 301         {
 302             settime[i] = readtime[i];
 303         }
 304     }
 305
 306     if (bSort && (nfiles > 1))
 307     {
 308         sort_files(fnms, settime, nfiles);
 309     }
 310     else
 311     {
 312         fprintf(stderr, "Sorting disabled.\n");
 313     }
 314
 315
 316     /* Write out the new order and start times */
 317     fprintf(stderr, "\nSummary of files and start times used:\n\n"
 318             "          File                Start time\n"
 319             "-----------------------------------------\n");
 320     for (i = 0; i < nfiles; i++)
 321     {
 322         switch (cont_type[i])
 323         {
 324             case TIME_EXPLICIT:
 325                 fprintf(stderr, "%25s   %10.3f\n", fnms[i], settime[i]);
 326                 break;
 327             case TIME_CONTINUE:
 328                 fprintf(stderr, "%25s        Continue from end of last file\n", fnms[i]);
 329                 break;
 330             case TIME_LAST:
 331                 fprintf(stderr, "%25s        Change by same amount as last file\n", fnms[i]);
 332                 break;
 333         }
 334     }
 335     fprintf(stderr, "\n");
 336
 337     settime[nfiles]   = FLT_MAX;
 338     cont_type[nfiles] = TIME_EXPLICIT;
 339     readtime[nfiles]  = FLT_MAX;
 340 }
 341
 342
 343 static void update_ee_sum(int nre,
 344                           gmx_int64_t *ee_sum_step,
 345                           gmx_int64_t *ee_sum_nsteps,
 346                           gmx_int64_t *ee_sum_nsum,
 347                           t_energy *ee_sum,
 348                           t_enxframe *fr, int out_step)
 349 {
 350     gmx_int64_t     nsteps, nsum, fr_nsum;
 351     int             i;
 352
 353     nsteps = *ee_sum_nsteps;
 354     nsum   = *ee_sum_nsum;
 355
 356     fr_nsum = fr->nsum;
 357     if (fr_nsum == 0)
 358     {
 359         fr_nsum = 1;
 360     }
 361
 362     if (nsteps == 0)
 363     {
 364         if (fr_nsum == 1)
 365         {
 366             for (i = 0; i < nre; i++)
 367             {
 368                 ee_sum[i].esum = fr->ener[i].e;
 369                 ee_sum[i].eav  = 0;
 370             }
 371         }
 372         else
 373         {
 374             for (i = 0; i < nre; i++)
 375             {
 376                 ee_sum[i].esum = fr->ener[i].esum;
 377                 ee_sum[i].eav  = fr->ener[i].eav;
 378             }
 379         }
 380         nsteps = fr->nsteps;
 381         nsum   = fr_nsum;
 382     }
 383     else if (out_step + *ee_sum_nsum - *ee_sum_step == nsteps + fr->nsteps)
 384     {
 385         if (fr_nsum == 1)
 386         {
 387             for (i = 0; i < nre; i++)
 388             {
 389                 ee_sum[i].eav  +=
 390                     gmx::square(ee_sum[i].esum/nsum
 391                                 - (ee_sum[i].esum + fr->ener[i].e)/(nsum + 1))*nsum*(nsum + 1);
 392                 ee_sum[i].esum += fr->ener[i].e;
 393             }
 394         }
 395         else
 396         {
 397             for (i = 0; i < fr->nre; i++)
 398             {
 399                 ee_sum[i].eav  +=
 400                     fr->ener[i].eav +
 401                     gmx::square(ee_sum[i].esum/nsum
 402                                 - (ee_sum[i].esum + fr->ener[i].esum)/(nsum + fr->nsum))*
 403                     nsum*(nsum + fr->nsum)/static_cast<double>(fr->nsum);
 404                 ee_sum[i].esum += fr->ener[i].esum;
 405             }
 406         }
 407         nsteps += fr->nsteps;
 408         nsum   += fr_nsum;
 409     }
 410     else
 411     {
 412         if (fr->nsum != 0)
 413         {
 414             fprintf(stderr, "\nWARNING: missing energy sums at time %f\n", fr->t);
 415         }
 416         nsteps = 0;
 417         nsum   = 0;
 418     }
 419
 420     *ee_sum_step   = out_step;
 421     *ee_sum_nsteps = nsteps;
 422     *ee_sum_nsum   = nsum;
 423 }
 424
 425 int gmx_eneconv(int argc, char *argv[])
 426 {
 427     const char       *desc[] = {
 428         "With [IT]multiple files[it] specified for the [TT]-f[tt] option:[PAR]",
 429         "Concatenates several energy files in sorted order.",
 430         "In the case of double time frames, the one",
 431         "in the later file is used. By specifying [TT]-settime[tt] you will be",
 432         "asked for the start time of each file. The input files are taken",
 433         "from the command line,",
 434         "such that the command [TT]gmx eneconv -f *.edr -o fixed.edr[tt] should do",
 435         "the trick. [PAR]",
 436         "With [IT]one file[it] specified for [TT]-f[tt]:[PAR]",
 437         "Reads one energy file and writes another, applying the [TT]-dt[tt],",
 438         "[TT]-offset[tt], [TT]-t0[tt] and [TT]-settime[tt] options and",
 439         "converting to a different format if necessary (indicated by file",
 440         "extentions).[PAR]",
 441         "[TT]-settime[tt] is applied first, then [TT]-dt[tt]/[TT]-offset[tt]",
 442         "followed by [TT]-b[tt] and [TT]-e[tt] to select which frames to write."
 443     };
 444     const char       *bugs[] = {
 445         "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."
 446     };
 447     ener_file_t       in  = NULL, out = NULL;
 448     gmx_enxnm_t      *enm = NULL;
 449 #if 0
 450     ener_file_t       in, out = NULL;
 451     gmx_enxnm_t      *enm = NULL;
 452 #endif
 453     t_enxframe       *fr, *fro;
 454     gmx_int64_t       ee_sum_step = 0, ee_sum_nsteps, ee_sum_nsum;
 455     t_energy         *ee_sum;
 456     gmx_int64_t       lastfilestep, laststep, startstep_file = 0;
 457     int               noutfr;
 458     int               nre, nremax, this_nre, nfile, f, i, kkk, nset, *set = NULL;
 459     double            last_t;
 460     char            **fnms;
 461     real             *readtime, *settime, timestep, tadjust;
 462     char              buf[22], buf2[22];
 463     int              *cont_type;
 464     gmx_bool          bNewFile, bFirst, bNewOutput;
 465     gmx_output_env_t *oenv;
 466     gmx_bool          warned_about_dh = FALSE;
 467     t_enxblock       *blocks          = NULL;
 468     int               nblocks         = 0;
 469     int               nblocks_alloc   = 0;
 470
 471     t_filenm          fnm[] = {
 472         { efEDR, "-f", NULL,    ffRDMULT },
 473         { efEDR, "-o", "fixed", ffWRITE  },
 474     };
 475
 476 #define NFILE asize(fnm)
 477     gmx_bool         bWrite;
 478     static real      delta_t   = 0.0, toffset = 0, scalefac = 1;
 479     static gmx_bool  bSetTime  = FALSE;
 480     static gmx_bool  bSort     = TRUE, bError = TRUE;
 481     static real      begin     = -1;
 482     static real      end       = -1;
 483     gmx_bool         remove_dh = FALSE;
 484
 485     t_pargs          pa[] = {
 486         { "-b",        FALSE, etREAL, {&begin},
 487           "First time to use"},
 488         { "-e",        FALSE, etREAL, {&end},
 489           "Last time to use"},
 490         { "-dt",       FALSE, etREAL, {&delta_t},
 491           "Only write out frame when t MOD dt = offset" },
 492         { "-offset",   FALSE, etREAL, {&toffset},
 493           "Time offset for [TT]-dt[tt] option" },
 494         { "-settime",  FALSE, etBOOL, {&bSetTime},
 495           "Change starting time interactively" },
 496         { "-sort",     FALSE, etBOOL, {&bSort},
 497           "Sort energy files (not frames)"},
 498         { "-rmdh",     FALSE, etBOOL, {&remove_dh},
 499           "Remove free energy block data" },
 500         { "-scalefac", FALSE, etREAL, {&scalefac},
 501           "Multiply energy component by this factor" },
 502         { "-error",    FALSE, etBOOL, {&bError},
 503           "Stop on errors in the file" }
 504     };
 505
 506     if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa),
 507                            pa, asize(desc), desc, asize(bugs), bugs, &oenv))
 508     {
 509         return 0;
 510     }
 511     tadjust  = 0;
 512     nremax   = 0;
 513     nset     = 0;
 514     timestep = 0.0;
 515     snew(fnms, argc);
 516     lastfilestep = 0;
 517     laststep     = 0;
 518
 519     nfile = opt2fns(&fnms, "-f", NFILE, fnm);
 520
 521     if (!nfile)
 522     {
 523         gmx_fatal(FARGS, "No input files!");
 524     }
 525
 526     snew(settime, nfile+1);
 527     snew(readtime, nfile+1);
 528     snew(cont_type, nfile+1);
 529
 530     nre = scan_ene_files(fnms, nfile, readtime, &timestep, &nremax);
 531     edit_files(fnms, nfile, 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 (f = 0; f < nfile; f++)
 548     {
 549         bNewFile   = TRUE;
 550         bNewOutput = TRUE;
 551         in         = open_enx(fnms[f], "r");
 552         enm        = NULL;
 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 %d\n",
 606                         gmx_step_str(fr->step, buf), fr->t,
 607                         gmx_step_str(fro->step, buf2), fro->t, bWrite);
 608             }
 609
 610             if (bError)
 611             {
 612                 if ((end > 0) && (fro->t > end+GMX_REAL_EPS))
 613                 {
 614                     f = nfile;
 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 = gmx_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                                                fnms[f], 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 == nfile)
 767         {
 768             f--;
 769         }
 770         printf("\nLast step written from %s: t %g, step %s\n",
 771                fnms[f], 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 < nfile-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 }