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