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