src/gromacs/gmxana/gmx_eneconv.c

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