src/gromacs/gmxana/gmx_eneconv.c

   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, 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 #ifdef HAVE_CONFIG_H
  38 #include <config.h>
  39 #endif
  40 #include <string.h>
  41 #include <math.h>
  42
  43 #include "typedefs.h"
  44 #include "gromacs/utility/smalloc.h"
  45 #include "gromacs/commandline/pargs.h"
  46 #include "disre.h"
  47 #include "names.h"
  48 #include "macros.h"
  49 #include "gmx_fatal.h"
  50 #include "gromacs/fileio/enxio.h"
  51 #include "vec.h"
  52 #include "gmx_ana.h"
  53 #include "gromacs/fileio/trxio.h"
  54
  55 #define TIME_EXPLICIT 0
  56 #define TIME_CONTINUE 1
  57 #define TIME_LAST     2
  58 #ifndef FLT_MAX
  59 #define FLT_MAX 1e36
  60 #endif
  61
  62 static int *select_it(int nre, gmx_enxnm_t *nm, int *nset)
  63 {
  64     gmx_bool *bE;
  65     int       n, k, j, i;
  66     int      *set;
  67     gmx_bool  bVerbose = TRUE;
  68
  69     if ((getenv("VERBOSE")) != NULL)
  70     {
  71         bVerbose = FALSE;
  72     }
  73
  74     fprintf(stderr, "Select the terms you want to scale from the following list\n");
  75     fprintf(stderr, "End your selection with 0\n");
  76
  77     if (bVerbose)
  78     {
  79         for (k = 0; (k < nre); )
  80         {
  81             for (j = 0; (j < 4) && (k < nre); j++, k++)
  82             {
  83                 fprintf(stderr, " %3d=%14s", k+1, nm[k].name);
  84             }
  85             fprintf(stderr, "\n");
  86         }
  87     }
  88
  89     snew(bE, nre);
  90     do
  91     {
  92         if (1 != scanf("%d", &n))
  93         {
  94             gmx_fatal(FARGS, "Cannot read energy term");
  95         }
  96         if ((n > 0) && (n <= nre))
  97         {
  98             bE[n-1] = TRUE;
  99         }
 100     }
 101     while (n != 0);
 102
 103     snew(set, nre);
 104     for (i = (*nset) = 0; (i < nre); i++)
 105     {
 106         if (bE[i])
 107         {
 108             set[(*nset)++] = i;
 109         }
 110     }
 111
 112     sfree(bE);
 113
 114     return set;
 115 }
 116
 117 static void sort_files(char **fnms, real *settime, int nfile)
 118 {
 119     int   i, j, minidx;
 120     real  timeswap;
 121     char *chptr;
 122
 123     for (i = 0; i < nfile; i++)
 124     {
 125         minidx = i;
 126         for (j = i+1; j < nfile; j++)
 127         {
 128             if (settime[j] < settime[minidx])
 129             {
 130                 minidx = j;
 131             }
 132         }
 133         if (minidx != i)
 134         {
 135             timeswap        = settime[i];
 136             settime[i]      = settime[minidx];
 137             settime[minidx] = timeswap;
 138             chptr           = fnms[i];
 139             fnms[i]         = fnms[minidx];
 140             fnms[minidx]    = chptr;
 141         }
 142     }
 143 }
 144
 145
 146 static int scan_ene_files(char **fnms, int nfiles,
 147                           real *readtime, real *timestep, int *nremax)
 148 {
 149     /* Check number of energy terms and start time of all files */
 150     int          f, i, nre, nremin = 0, nresav = 0;
 151     ener_file_t  in;
 152     real         t1, t2;
 153     char         inputstring[STRLEN];
 154     gmx_enxnm_t *enm;
 155     t_enxframe  *fr;
 156
 157     snew(fr, 1);
 158
 159     for (f = 0; f < nfiles; f++)
 160     {
 161         in  = open_enx(fnms[f], "r");
 162         enm = NULL;
 163         do_enxnms(in, &nre, &enm);
 164
 165         if (f == 0)
 166         {
 167             nresav  = nre;
 168             nremin  = nre;
 169             *nremax = nre;
 170             do_enx(in, fr);
 171             t1 = fr->t;
 172             do_enx(in, fr);
 173             t2          = fr->t;
 174             *timestep   = t2-t1;
 175             readtime[f] = t1;
 176             close_enx(in);
 177         }
 178         else
 179         {
 180             nremin  = min(nremin, fr->nre);
 181             *nremax = max(*nremax, fr->nre);
 182             if (nre != nresav)
 183             {
 184                 fprintf(stderr,
 185                         "Energy files don't match, different number of energies:\n"
 186                         " %s: %d\n %s: %d\n", fnms[f-1], nresav, fnms[f], fr->nre);
 187                 fprintf(stderr,
 188                         "\nContinue conversion using only the first %d terms (n/y)?\n"
 189                         "(you should be sure that the energy terms match)\n", nremin);
 190                 if (NULL == fgets(inputstring, STRLEN-1, stdin))
 191                 {
 192                     gmx_fatal(FARGS, "Error reading user input");
 193                 }
 194                 if (inputstring[0] != 'y' && inputstring[0] != 'Y')
 195                 {
 196                     fprintf(stderr, "Will not convert\n");
 197                     exit(0);
 198                 }
 199                 nresav = fr->nre;
 200             }
 201             do_enx(in, fr);
 202             readtime[f] = fr->t;
 203             close_enx(in);
 204         }
 205         fprintf(stderr, "\n");
 206         free_enxnms(nre, enm);
 207     }
 208
 209     free_enxframe(fr);
 210     sfree(fr);
 211
 212     return nremin;
 213 }
 214
 215
 216 static void edit_files(char **fnms, int nfiles, real *readtime,
 217                        real *settime, int *cont_type, gmx_bool bSetTime, gmx_bool bSort)
 218 {
 219     int      i;
 220     gmx_bool ok;
 221     char     inputstring[STRLEN], *chptr;
 222
 223     if (bSetTime)
 224     {
 225         if (nfiles == 1)
 226         {
 227             fprintf(stderr, "\n\nEnter the new start time:\n\n");
 228         }
 229         else
 230         {
 231             fprintf(stderr, "\n\nEnter the new start time for each file.\n"
 232                     "There are two special options, both disables sorting:\n\n"
 233                     "c (continue) - The start time is taken from the end\n"
 234                     "of the previous file. Use it when your continuation run\n"
 235                     "restarts with t=0 and there is no overlap.\n\n"
 236                     "l (last) - The time in this file will be changed the\n"
 237                     "same amount as in the previous. Use it when the time in the\n"
 238                     "new run continues from the end of the previous one,\n"
 239                     "since this takes possible overlap into account.\n\n");
 240         }
 241
 242         fprintf(stderr, "          File             Current start       New start\n"
 243                 "---------------------------------------------------------\n");
 244
 245         for (i = 0; i < nfiles; i++)
 246         {
 247             fprintf(stderr, "%25s   %10.3f             ", fnms[i], readtime[i]);
 248             ok = FALSE;
 249             do
 250             {
 251                 if (NULL == fgets(inputstring, STRLEN-1, stdin))
 252                 {
 253                     gmx_fatal(FARGS, "Error reading user input");
 254                 }
 255                 inputstring[strlen(inputstring)-1] = 0;
 256
 257                 if (inputstring[0] == 'c' || inputstring[0] == 'C')
 258                 {
 259                     cont_type[i] = TIME_CONTINUE;
 260                     bSort        = FALSE;
 261                     ok           = TRUE;
 262                     settime[i]   = FLT_MAX;
 263                 }
 264                 else if (inputstring[0] == 'l' ||
 265                          inputstring[0] == 'L')
 266                 {
 267                     cont_type[i] = TIME_LAST;
 268                     bSort        = FALSE;
 269                     ok           = TRUE;
 270                     settime[i]   = FLT_MAX;
 271                 }
 272                 else
 273                 {
 274                     settime[i] = strtod(inputstring, &chptr);
 275                     if (chptr == inputstring)
 276                     {
 277                         fprintf(stderr, "Try that again: ");
 278                     }
 279                     else
 280                     {
 281                         cont_type[i] = TIME_EXPLICIT;
 282                         ok           = TRUE;
 283                     }
 284                 }
 285             }
 286             while (!ok);
 287         }
 288         if (cont_type[0] != TIME_EXPLICIT)
 289         {
 290             cont_type[0] = TIME_EXPLICIT;
 291             settime[0]   = 0;
 292         }
 293     }
 294     else
 295     {
 296         for (i = 0; i < nfiles; i++)
 297         {
 298             settime[i] = readtime[i];
 299         }
 300     }
 301
 302     if (bSort && (nfiles > 1))
 303     {
 304         sort_files(fnms, settime, nfiles);
 305     }
 306     else
 307     {
 308         fprintf(stderr, "Sorting disabled.\n");
 309     }
 310
 311
 312     /* Write out the new order and start times */
 313     fprintf(stderr, "\nSummary of files and start times used:\n\n"
 314             "          File                Start time\n"
 315             "-----------------------------------------\n");
 316     for (i = 0; i < nfiles; i++)
 317     {
 318         switch (cont_type[i])
 319         {
 320             case TIME_EXPLICIT:
 321                 fprintf(stderr, "%25s   %10.3f\n", fnms[i], settime[i]);
 322                 break;
 323             case TIME_CONTINUE:
 324                 fprintf(stderr, "%25s        Continue from end of last file\n", fnms[i]);
 325                 break;
 326             case TIME_LAST:
 327                 fprintf(stderr, "%25s        Change by same amount as last file\n", fnms[i]);
 328                 break;
 329         }
 330     }
 331     fprintf(stderr, "\n");
 332
 333     settime[nfiles]   = FLT_MAX;
 334     cont_type[nfiles] = TIME_EXPLICIT;
 335     readtime[nfiles]  = FLT_MAX;
 336 }
 337
 338
 339 static void copy_ee(t_energy *src, t_energy *dst, int nre)
 340 {
 341     int i;
 342
 343     for (i = 0; i < nre; i++)
 344     {
 345         dst[i].e    = src[i].e;
 346         dst[i].esum = src[i].esum;
 347         dst[i].eav  = src[i].eav;
 348     }
 349 }
 350
 351 static void update_ee(t_energy *lastee, gmx_int64_t laststep,
 352                       t_energy *startee, gmx_int64_t startstep,
 353                       t_energy *ee, int step,
 354                       t_energy *outee, int nre)
 355 {
 356     int    i;
 357     double sigmacorr, nom, denom;
 358     double prestart_esum;
 359     double prestart_sigma;
 360
 361     for (i = 0; i < nre; i++)
 362     {
 363         outee[i].e = ee[i].e;
 364         /* add statistics from earlier file if present */
 365         if (laststep > 0)
 366         {
 367             outee[i].esum = lastee[i].esum+ee[i].esum;
 368             nom           = (lastee[i].esum*(step+1)-ee[i].esum*(laststep));
 369             denom         = ((step+1.0)*(laststep)*(step+1.0+laststep));
 370             sigmacorr     = nom*nom/denom;
 371             outee[i].eav  = lastee[i].eav+ee[i].eav+sigmacorr;
 372         }
 373         else
 374         {
 375             /* otherwise just copy to output */
 376             outee[i].esum = ee[i].esum;
 377             outee[i].eav  = ee[i].eav;
 378         }
 379
 380         /* if we didnt start to write at the first frame
 381          * we must compensate the statistics for this
 382          * there are laststep frames in the earlier file
 383          * and step+1 frames in this one.
 384          */
 385         if (startstep > 0)
 386         {
 387             gmx_int64_t q = laststep+step;
 388             gmx_int64_t p = startstep+1;
 389             prestart_esum  = startee[i].esum-startee[i].e;
 390             sigmacorr      = prestart_esum-(p-1)*startee[i].e;
 391             prestart_sigma = startee[i].eav-
 392                 sigmacorr*sigmacorr/(p*(p-1));
 393             sigmacorr = prestart_esum/(p-1)-
 394                 outee[i].esum/(q);
 395             outee[i].esum -= prestart_esum;
 396             if (q-p+1 > 0)
 397             {
 398                 outee[i].eav = outee[i].eav-prestart_sigma-
 399                     sigmacorr*sigmacorr*((p-1)*q)/(q-p+1);
 400             }
 401         }
 402
 403         if ((outee[i].eav/(laststep+step+1)) < (GMX_REAL_EPS))
 404         {
 405             outee[i].eav = 0;
 406         }
 407     }
 408 }
 409
 410 static void update_ee_sum(int nre,
 411                           gmx_int64_t *ee_sum_step,
 412                           gmx_int64_t *ee_sum_nsteps,
 413                           gmx_int64_t *ee_sum_nsum,
 414                           t_energy *ee_sum,
 415                           t_enxframe *fr, int out_step)
 416 {
 417     gmx_int64_t     nsteps, nsum, fr_nsum;
 418     int             i;
 419
 420     nsteps = *ee_sum_nsteps;
 421     nsum   = *ee_sum_nsum;
 422
 423     fr_nsum = fr->nsum;
 424     if (fr_nsum == 0)
 425     {
 426         fr_nsum = 1;
 427     }
 428
 429     if (nsteps == 0)
 430     {
 431         if (fr_nsum == 1)
 432         {
 433             for (i = 0; i < nre; i++)
 434             {
 435                 ee_sum[i].esum = fr->ener[i].e;
 436                 ee_sum[i].eav  = 0;
 437             }
 438         }
 439         else
 440         {
 441             for (i = 0; i < nre; i++)
 442             {
 443                 ee_sum[i].esum = fr->ener[i].esum;
 444                 ee_sum[i].eav  = fr->ener[i].eav;
 445             }
 446         }
 447         nsteps = fr->nsteps;
 448         nsum   = fr_nsum;
 449     }
 450     else if (out_step + *ee_sum_nsum - *ee_sum_step == nsteps + fr->nsteps)
 451     {
 452         if (fr_nsum == 1)
 453         {
 454             for (i = 0; i < nre; i++)
 455             {
 456                 ee_sum[i].eav  +=
 457                     dsqr(ee_sum[i].esum/nsum
 458                          - (ee_sum[i].esum + fr->ener[i].e)/(nsum + 1))*nsum*(nsum + 1);
 459                 ee_sum[i].esum += fr->ener[i].e;
 460             }
 461         }
 462         else
 463         {
 464             for (i = 0; i < fr->nre; i++)
 465             {
 466                 ee_sum[i].eav  +=
 467                     fr->ener[i].eav +
 468                     dsqr(ee_sum[i].esum/nsum
 469                          - (ee_sum[i].esum + fr->ener[i].esum)/(nsum + fr->nsum))*
 470                     nsum*(nsum + fr->nsum)/(double)fr->nsum;
 471                 ee_sum[i].esum += fr->ener[i].esum;
 472             }
 473         }
 474         nsteps += fr->nsteps;
 475         nsum   += fr_nsum;
 476     }
 477     else
 478     {
 479         if (fr->nsum != 0)
 480         {
 481             fprintf(stderr, "\nWARNING: missing energy sums at time %f\n", fr->t);
 482         }
 483         nsteps = 0;
 484         nsum   = 0;
 485     }
 486
 487     *ee_sum_step   = out_step;
 488     *ee_sum_nsteps = nsteps;
 489     *ee_sum_nsum   = nsum;
 490 }
 491
 492 int gmx_eneconv(int argc, char *argv[])
 493 {
 494     const char     *desc[] = {
 495         "With [IT]multiple files[it] specified for the [TT]-f[tt] option:[BR]",
 496         "Concatenates several energy files in sorted order.",
 497         "In the case of double time frames, the one",
 498         "in the later file is used. By specifying [TT]-settime[tt] you will be",
 499         "asked for the start time of each file. The input files are taken",
 500         "from the command line,",
 501         "such that the command [TT]gmx eneconv -f *.edr -o fixed.edr[tt] should do",
 502         "the trick. [PAR]",
 503         "With [IT]one file[it] specified for [TT]-f[tt]:[BR]",
 504         "Reads one energy file and writes another, applying the [TT]-dt[tt],",
 505         "[TT]-offset[tt], [TT]-t0[tt] and [TT]-settime[tt] options and",
 506         "converting to a different format if necessary (indicated by file",
 507         "extentions).[PAR]",
 508         "[TT]-settime[tt] is applied first, then [TT]-dt[tt]/[TT]-offset[tt]",
 509         "followed by [TT]-b[tt] and [TT]-e[tt] to select which frames to write."
 510     };
 511     const char     *bugs[] = {
 512         "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."
 513     };
 514     ener_file_t     in  = NULL, out = NULL;
 515     gmx_enxnm_t    *enm = NULL;
 516 #if 0
 517     ener_file_t     in, out = NULL;
 518     gmx_enxnm_t    *enm = NULL;
 519 #endif
 520     t_enxframe     *fr, *fro;
 521     gmx_int64_t     ee_sum_step = 0, ee_sum_nsteps, ee_sum_nsum;
 522     t_energy       *ee_sum;
 523     gmx_int64_t     lastfilestep, laststep, startstep, startstep_file = 0;
 524     int             noutfr;
 525     int             nre, nremax, this_nre, nfile, f, i, j, kkk, nset, *set = NULL;
 526     double          last_t;
 527     char          **fnms;
 528     real           *readtime, *settime, timestep, t1, tadjust;
 529     char            inputstring[STRLEN], *chptr, buf[22], buf2[22], buf3[22];
 530     gmx_bool        ok;
 531     int            *cont_type;
 532     gmx_bool        bNewFile, bFirst, bNewOutput;
 533     output_env_t    oenv;
 534     gmx_bool        warned_about_dh = FALSE;
 535     t_enxblock     *blocks          = NULL;
 536     int             nblocks         = 0;
 537     int             nblocks_alloc   = 0;
 538
 539     t_filenm        fnm[] = {
 540         { efEDR, "-f", NULL,    ffRDMULT },
 541         { efEDR, "-o", "fixed", ffWRITE  },
 542     };
 543
 544 #define NFILE asize(fnm)
 545     gmx_bool         bWrite;
 546     static real      delta_t   = 0.0, toffset = 0, scalefac = 1;
 547     static gmx_bool  bSetTime  = FALSE;
 548     static gmx_bool  bSort     = TRUE, bError = TRUE;
 549     static real      begin     = -1;
 550     static real      end       = -1;
 551     gmx_bool         remove_dh = FALSE;
 552
 553     t_pargs          pa[] = {
 554         { "-b",        FALSE, etREAL, {&begin},
 555           "First time to use"},
 556         { "-e",        FALSE, etREAL, {&end},
 557           "Last time to use"},
 558         { "-dt",       FALSE, etREAL, {&delta_t},
 559           "Only write out frame when t MOD dt = offset" },
 560         { "-offset",   FALSE, etREAL, {&toffset},
 561           "Time offset for [TT]-dt[tt] option" },
 562         { "-settime",  FALSE, etBOOL, {&bSetTime},
 563           "Change starting time interactively" },
 564         { "-sort",     FALSE, etBOOL, {&bSort},
 565           "Sort energy files (not frames)"},
 566         { "-rmdh",     FALSE, etBOOL, {&remove_dh},
 567           "Remove free energy block data" },
 568         { "-scalefac", FALSE, etREAL, {&scalefac},
 569           "Multiply energy component by this factor" },
 570         { "-error",    FALSE, etBOOL, {&bError},
 571           "Stop on errors in the file" }
 572     };
 573
 574     if (!parse_common_args(&argc, argv, PCA_BE_NICE, NFILE, fnm, asize(pa),
 575                            pa, asize(desc), desc, asize(bugs), bugs, &oenv))
 576     {
 577         return 0;
 578     }
 579     tadjust  = 0;
 580     nremax   = 0;
 581     nset     = 0;
 582     timestep = 0.0;
 583     snew(fnms, argc);
 584     nfile        = 0;
 585     lastfilestep = 0;
 586     laststep     = startstep = 0;
 587
 588     nfile = opt2fns(&fnms, "-f", NFILE, fnm);
 589
 590     if (!nfile)
 591     {
 592         gmx_fatal(FARGS, "No input files!");
 593     }
 594
 595     snew(settime, nfile+1);
 596     snew(readtime, nfile+1);
 597     snew(cont_type, nfile+1);
 598
 599     nre = scan_ene_files(fnms, nfile, readtime, &timestep, &nremax);
 600     edit_files(fnms, nfile, readtime, settime, cont_type, bSetTime, bSort);
 601
 602     ee_sum_nsteps = 0;
 603     ee_sum_nsum   = 0;
 604     snew(ee_sum, nremax);
 605
 606     snew(fr, 1);
 607     snew(fro, 1);
 608     fro->t   = -1e20;
 609     fro->nre = nre;
 610     snew(fro->ener, nremax);
 611
 612     noutfr = 0;
 613     bFirst = TRUE;
 614
 615     last_t = fro->t;
 616     for (f = 0; f < nfile; f++)
 617     {
 618         bNewFile   = TRUE;
 619         bNewOutput = TRUE;
 620         in         = open_enx(fnms[f], "r");
 621         enm        = NULL;
 622         do_enxnms(in, &this_nre, &enm);
 623         if (f == 0)
 624         {
 625             if (scalefac != 1)
 626             {
 627                 set = select_it(nre, enm, &nset);
 628             }
 629
 630             /* write names to the output file */
 631             out = open_enx(opt2fn("-o", NFILE, fnm), "w");
 632             do_enxnms(out, &nre, &enm);
 633         }
 634
 635         /* start reading from the next file */
 636         while ((fro->t <= (settime[f+1] + GMX_REAL_EPS)) &&
 637                do_enx(in, fr))
 638         {
 639             if (bNewFile)
 640             {
 641                 startstep_file = fr->step;
 642                 tadjust        = settime[f] - fr->t;
 643                 if (cont_type[f+1] == TIME_LAST)
 644                 {
 645                     settime[f+1]   = readtime[f+1]-readtime[f]+settime[f];
 646                     cont_type[f+1] = TIME_EXPLICIT;
 647                 }
 648                 bNewFile = FALSE;
 649             }
 650
 651             if (tadjust + fr->t <= last_t)
 652             {
 653                 /* Skip this frame, since we already have it / past it */
 654                 if (debug)
 655                 {
 656                     fprintf(debug, "fr->step %s, fr->t %.4f\n",
 657                             gmx_step_str(fr->step, buf), fr->t);
 658                     fprintf(debug, "tadjust %12.6e + fr->t %12.6e <= t %12.6e\n",
 659                             tadjust, fr->t, last_t);
 660                 }
 661                 continue;
 662             }
 663
 664             fro->step = lastfilestep + fr->step - startstep_file;
 665             fro->t    = tadjust  + fr->t;
 666
 667             bWrite = ((begin < 0 || (begin >= 0 && (fro->t >= begin-GMX_REAL_EPS))) &&
 668                       (end  < 0 || (end  >= 0 && (fro->t <= end  +GMX_REAL_EPS))) &&
 669                       (fro->t <= settime[f+1]+0.5*timestep));
 670
 671             if (debug)
 672             {
 673                 fprintf(debug,
 674                         "fr->step %s, fr->t %.4f, fro->step %s fro->t %.4f, w %d\n",
 675                         gmx_step_str(fr->step, buf), fr->t,
 676                         gmx_step_str(fro->step, buf2), fro->t, bWrite);
 677             }
 678
 679             if (bError)
 680             {
 681                 if ((end > 0) && (fro->t > end+GMX_REAL_EPS))
 682                 {
 683                     f = nfile;
 684                     break;
 685                 }
 686             }
 687
 688             if (fro->t >= begin-GMX_REAL_EPS)
 689             {
 690                 if (bFirst)
 691                 {
 692                     bFirst    = FALSE;
 693                     startstep = fr->step;
 694                 }
 695                 if (bWrite)
 696                 {
 697                     update_ee_sum(nre, &ee_sum_step, &ee_sum_nsteps, &ee_sum_nsum, ee_sum,
 698                                   fr, fro->step);
 699                 }
 700             }
 701
 702             /* determine if we should write it */
 703             if (bWrite && (delta_t == 0 || bRmod(fro->t, toffset, delta_t)))
 704             {
 705                 laststep = fro->step;
 706                 last_t   = fro->t;
 707                 if (bNewOutput)
 708                 {
 709                     bNewOutput = FALSE;
 710                     fprintf(stderr, "\nContinue writing frames from t=%g, step=%s\n",
 711                             fro->t, gmx_step_str(fro->step, buf));
 712                 }
 713
 714                 /* Copy the energies */
 715                 for (i = 0; i < nre; i++)
 716                 {
 717                     fro->ener[i].e = fr->ener[i].e;
 718                 }
 719
 720                 fro->nsteps = ee_sum_nsteps;
 721                 fro->dt     = fr->dt;
 722
 723                 if (ee_sum_nsum <= 1)
 724                 {
 725                     fro->nsum = 0;
 726                 }
 727                 else
 728                 {
 729                     fro->nsum = gmx_int64_to_int(ee_sum_nsum,
 730                                                  "energy average summation");
 731                     /* Copy the energy sums */
 732                     for (i = 0; i < nre; i++)
 733                     {
 734                         fro->ener[i].esum = ee_sum[i].esum;
 735                         fro->ener[i].eav  = ee_sum[i].eav;
 736                     }
 737                 }
 738                 /* We wrote the energies, so reset the counts */
 739                 ee_sum_nsteps = 0;
 740                 ee_sum_nsum   = 0;
 741
 742                 if (scalefac != 1)
 743                 {
 744                     for (kkk = 0; kkk < nset; kkk++)
 745                     {
 746                         fro->ener[set[kkk]].e    *= scalefac;
 747                         if (fro->nsum > 0)
 748                         {
 749                             fro->ener[set[kkk]].eav  *= scalefac*scalefac;
 750                             fro->ener[set[kkk]].esum *= scalefac;
 751                         }
 752                     }
 753                 }
 754                 /* Copy restraint stuff */
 755                 /*fro->ndisre       = fr->ndisre;
 756                    fro->disre_rm3tav = fr->disre_rm3tav;
 757                    fro->disre_rt     = fr->disre_rt;*/
 758                 fro->nblock       = fr->nblock;
 759                 /*fro->nr           = fr->nr;*/
 760                 fro->block        = fr->block;
 761
 762                 /* check if we have blocks with delta_h data and are throwing
 763                    away data */
 764                 if (fro->nblock > 0)
 765                 {
 766                     if (remove_dh)
 767                     {
 768                         int i;
 769                         if (!blocks || nblocks_alloc < fr->nblock)
 770                         {
 771                             /* we pre-allocate the blocks */
 772                             nblocks_alloc = fr->nblock;
 773                             snew(blocks, nblocks_alloc);
 774                         }
 775                         nblocks = 0; /* number of blocks so far */
 776
 777                         for (i = 0; i < fr->nblock; i++)
 778                         {
 779                             if ( (fr->block[i].id != enxDHCOLL) &&
 780                                  (fr->block[i].id != enxDH) &&
 781                                  (fr->block[i].id != enxDHHIST) )
 782                             {
 783                                 /* copy everything verbatim */
 784                                 blocks[nblocks] = fr->block[i];
 785                                 nblocks++;
 786                             }
 787                         }
 788                         /* now set the block pointer to the new blocks */
 789                         fro->nblock = nblocks;
 790                         fro->block  = blocks;
 791                     }
 792                     else if (delta_t > 0)
 793                     {
 794                         if (!warned_about_dh)
 795                         {
 796                             for (i = 0; i < fr->nblock; i++)
 797                             {
 798                                 if (fr->block[i].id == enxDH ||
 799                                     fr->block[i].id == enxDHHIST)
 800                                 {
 801                                     int size;
 802                                     if (fr->block[i].id == enxDH)
 803                                     {
 804                                         size = fr->block[i].sub[2].nr;
 805                                     }
 806                                     else
 807                                     {
 808                                         size = fr->nsteps;
 809                                     }
 810                                     if (size > 0)
 811                                     {
 812                                         printf("\nWARNING: %s contains delta H blocks or histograms for which\n"
 813                                                "         some data is thrown away on a block-by-block basis, where each block\n"
 814                                                "         contains up to %d samples.\n"
 815                                                "         This is almost certainly not what you want.\n"
 816                                                "         Use the -rmdh option to throw all delta H samples away.\n"
 817                                                "         Use g_energy -odh option to extract these samples.\n",
 818                                                fnms[f], size);
 819                                         warned_about_dh = TRUE;
 820                                         break;
 821                                     }
 822                                 }
 823                             }
 824                         }
 825                     }
 826                 }
 827
 828                 do_enx(out, fro);
 829                 if (noutfr % 1000 == 0)
 830                 {
 831                     fprintf(stderr, "Writing frame time %g    ", fro->t);
 832                 }
 833                 noutfr++;
 834             }
 835         }
 836         if (f == nfile)
 837         {
 838             f--;
 839         }
 840         printf("\nLast step written from %s: t %g, step %s\n",
 841                fnms[f], last_t, gmx_step_str(laststep, buf));
 842         lastfilestep = laststep;
 843
 844         /* set the next time from the last in previous file */
 845         if (cont_type[f+1] == TIME_CONTINUE)
 846         {
 847             settime[f+1] = fro->t;
 848             /* in this case we have already written the last frame of
 849              * previous file, so update begin to avoid doubling it
 850              * with the start of the next file
 851              */
 852             begin = fro->t+0.5*timestep;
 853             /* cont_type[f+1]==TIME_EXPLICIT; */
 854         }
 855
 856         if ((fro->t < end) && (f < nfile-1) &&
 857             (fro->t < settime[f+1]-1.5*timestep))
 858         {
 859             fprintf(stderr,
 860                     "\nWARNING: There might be a gap around t=%g\n", fro->t);
 861         }
 862
 863         /* move energies to lastee */
 864         close_enx(in);
 865         free_enxnms(this_nre, enm);
 866
 867         fprintf(stderr, "\n");
 868     }
 869     if (noutfr == 0)
 870     {
 871         fprintf(stderr, "No frames written.\n");
 872     }
 873     else
 874     {
 875         fprintf(stderr, "Last frame written was at step %s, time %f\n",
 876                 gmx_step_str(fro->step, buf), fro->t);
 877         fprintf(stderr, "Wrote %d frames\n", noutfr);
 878     }
 879
 880     return 0;
 881 }