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