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