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