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   bool *bE;
  64   int  n,k,j,i;
  65   int  *set;
  66   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 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 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,bool bSetTime,bool bSort)
 210 {
 211   int i;
 212   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 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 -o fixed.edr *.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   bool       ok;
 500   int        *cont_type;
 501   bool       bNewFile,bFirst,bNewOutput;
 502   output_env_t oenv;
 503
 504   t_filenm fnm[] = {
 505     { efEDR, "-f", NULL,    ffRDMULT },
 506     { efEDR, "-o", "fixed", ffWRITE  },
 507   };
 508
 509 #define NFILE asize(fnm)
 510   bool   bWrite;
 511   static real  delta_t=0.0, toffset=0,scalefac=1;
 512   static bool  bSetTime=FALSE;
 513   static bool  bSort=TRUE,bError=TRUE;
 514   static real  begin=-1;
 515   static real  end=-1;
 516
 517   t_pargs pa[] = {
 518     { "-b",        FALSE, etREAL, {&begin},
 519       "First time to use"},
 520     { "-e",        FALSE, etREAL, {&end},
 521       "Last time to use"},
 522     { "-dt",       FALSE, etREAL, {&delta_t},
 523       "Only write out frame when t MOD dt = offset" },
 524     { "-offset",   FALSE, etREAL, {&toffset},
 525       "Time offset for -dt option" },
 526     { "-settime",  FALSE, etBOOL, {&bSetTime},
 527       "Change starting time interactively" },
 528     { "-sort",     FALSE, etBOOL, {&bSort},
 529       "Sort energy files (not frames)"},
 530     { "-scalefac", FALSE, etREAL, {&scalefac},
 531       "Multiply energy component by this factor" },
 532     { "-error",    FALSE, etBOOL, {&bError},
 533       "Stop on errors in the file" }
 534   };
 535
 536   CopyRight(stderr,argv[0]);
 537   parse_common_args(&argc,argv,PCA_BE_NICE,NFILE,fnm,asize(pa),
 538                     pa,asize(desc),desc,asize(bugs),bugs,&oenv);
 539   tadjust  = 0;
 540   nremax   = 0;
 541   nset     = 0;
 542   timestep = 0.0;
 543   snew(fnms,argc);
 544   nfile=0;
 545   lastfilestep=0;
 546   laststep=startstep=0;
 547
 548   nfile = opt2fns(&fnms,"-f",NFILE,fnm);
 549
 550   if (!nfile)
 551     gmx_fatal(FARGS,"No input files!");
 552
 553   snew(settime,nfile+1);
 554   snew(readtime,nfile+1);
 555   snew(cont_type,nfile+1);
 556
 557   nre=scan_ene_files(fnms,nfile,readtime,&timestep,&nremax);
 558   edit_files(fnms,nfile,readtime,settime,cont_type,bSetTime,bSort);
 559
 560   ee_sum_nsteps = 0;
 561   ee_sum_nsum   = 0;
 562   snew(ee_sum,nremax);
 563
 564   snew(fr,1);
 565   snew(fro,1);
 566   fro->t = -1e20;
 567   fro->nre = nre;
 568   snew(fro->ener,nremax);
 569
 570   noutfr=0;
 571   bFirst=TRUE;
 572
 573   last_t = fro->t;
 574   for(f=0;f<nfile;f++) {
 575     bNewFile=TRUE;
 576     bNewOutput=TRUE;
 577     in=open_enx(fnms[f],"r");
 578     enm = NULL;
 579     do_enxnms(in,&this_nre,&enm);
 580     if (f == 0) {
 581       if (scalefac != 1)
 582         set = select_it(nre,enm,&nset);
 583
 584       /* write names to the output file */
 585       out=open_enx(opt2fn("-o",NFILE,fnm),"w");
 586       do_enxnms(out,&nre,&enm);
 587     }
 588
 589     /* start reading from the next file */
 590     while ((fro->t <= (settime[f+1] + GMX_REAL_EPS)) &&
 591           do_enx(in,fr)) {
 592       if(bNewFile) {
 593         startstep_file = fr->step;
 594         tadjust = settime[f] - fr->t;
 595         if(cont_type[f+1]==TIME_LAST) {
 596           settime[f+1]   = readtime[f+1]-readtime[f]+settime[f];
 597           cont_type[f+1] = TIME_EXPLICIT;
 598         }
 599         bNewFile = FALSE;
 600       }
 601
 602       if (tadjust + fr->t <= last_t) {
 603         /* Skip this frame, since we already have it / past it */
 604         if (debug) {
 605           fprintf(debug,"fr->step %s, fr->t %.4f\n",
 606                   gmx_step_str(fr->step,buf),fr->t);
 607           fprintf(debug,"tadjust %12.6e + fr->t %12.6e <= t %12.6e\n",
 608                   tadjust,fr->t,last_t);
 609         }
 610         continue;
 611       }
 612
 613       fro->step = lastfilestep + fr->step - startstep_file;
 614       fro->t    = tadjust  + fr->t;
 615
 616       bWrite = ((begin<0 || (begin>=0 && (fro->t >= begin-GMX_REAL_EPS))) &&
 617                 (end  <0 || (end  >=0 && (fro->t <= end  +GMX_REAL_EPS))) &&
 618                 (fro->t <= settime[f+1]+0.5*timestep));
 619
 620       if (debug) {
 621         fprintf(debug,
 622                 "fr->step %s, fr->t %.4f, fro->step %s fro->t %.4f, w %d\n",
 623                 gmx_step_str(fr->step,buf),fr->t,
 624                 gmx_step_str(fro->step,buf2),fro->t,bWrite);
 625       }
 626
 627       if (bError)
 628         if ((end > 0) && (fro->t > end+GMX_REAL_EPS)) {
 629           f = nfile;
 630           break;
 631         }
 632
 633       if (fro->t >= begin-GMX_REAL_EPS) {
 634         if (bFirst) {
 635           bFirst = FALSE;
 636           startstep = fr->step;
 637         }
 638         if (bWrite) {
 639           update_ee_sum(nre,&ee_sum_step,&ee_sum_nsteps,&ee_sum_nsum,ee_sum,
 640                         fr,fro->step);
 641         }
 642       }
 643
 644       /* determine if we should write it */
 645       if (bWrite && (delta_t==0 || bRmod(fro->t,toffset,delta_t))) {
 646         laststep = fro->step;
 647         last_t   = fro->t;
 648         if(bNewOutput) {
 649           bNewOutput=FALSE;
 650           fprintf(stderr,"\nContinue writing frames from t=%g, step=%s\n",
 651                   fro->t,gmx_step_str(fro->step,buf));
 652         }
 653
 654         /* Copy the energies */
 655         for(i=0; i<nre; i++) {
 656           fro->ener[i].e = fr->ener[i].e;
 657         }
 658
 659         fro->nsteps = ee_sum_nsteps;
 660         fro->dt     = fr->dt;
 661
 662         if (ee_sum_nsum <= 1) {
 663           fro->nsum = 0;
 664         } else {
 665           fro->nsum = gmx_large_int_to_int(ee_sum_nsum,
 666                                            "energy average summation");
 667           /* Copy the energy sums */
 668           for(i=0; i<nre; i++) {
 669             fro->ener[i].esum = ee_sum[i].esum;
 670             fro->ener[i].eav  = ee_sum[i].eav;
 671           }
 672         }
 673         /* We wrote the energies, so reset the counts */
 674         ee_sum_nsteps = 0;
 675         ee_sum_nsum   = 0;
 676
 677         if (scalefac != 1) {
 678           for(kkk=0; kkk<nset; kkk++) {
 679             fro->ener[set[kkk]].e    *= scalefac;
 680             if (fro->nsum > 0) {
 681               fro->ener[set[kkk]].eav  *= scalefac*scalefac;
 682               fro->ener[set[kkk]].esum *= scalefac;
 683             }
 684           }
 685         }
 686         /* Copy restraint stuff */
 687         /*fro->ndisre       = fr->ndisre;
 688         fro->disre_rm3tav = fr->disre_rm3tav;
 689         fro->disre_rt     = fr->disre_rt;*/
 690         fro->nblock       = fr->nblock;
 691         /*fro->nr           = fr->nr;*/
 692         fro->block        = fr->block;
 693
 694         do_enx(out,fro);
 695         if (noutfr % 1000 == 0)
 696           fprintf(stderr,"Writing frame time %g    ",fro->t);
 697         noutfr++;
 698       }
 699     }
 700
 701     printf("\nLast step written from %s: t %g, step %s\n",
 702            fnms[f],last_t,gmx_step_str(laststep,buf));
 703     lastfilestep = laststep;
 704
 705     /* set the next time from the last in previous file */
 706     if (cont_type[f+1]==TIME_CONTINUE) {
 707         settime[f+1] = fro->t;
 708         /* in this case we have already written the last frame of
 709          * previous file, so update begin to avoid doubling it
 710          * with the start of the next file
 711          */
 712         begin = fro->t+0.5*timestep;
 713         /* cont_type[f+1]==TIME_EXPLICIT; */
 714     }
 715
 716     if ((fro->t < end) && (f < nfile-1) &&
 717         (fro->t < settime[f+1]-1.5*timestep))
 718       fprintf(stderr,
 719               "\nWARNING: There might be a gap around t=%g\n",fro->t);
 720
 721     /* move energies to lastee */
 722     close_enx(in);
 723     free_enxnms(this_nre,enm);
 724
 725     fprintf(stderr,"\n");
 726   }
 727   if (noutfr == 0)
 728     fprintf(stderr,"No frames written.\n");
 729   else {
 730     fprintf(stderr,"Last frame written was at step %s, time %f\n",
 731             gmx_step_str(fro->step,buf),fro->t);
 732     fprintf(stderr,"Wrote %d frames\n",noutfr);
 733   }
 734
 735   thanx(stderr);
 736   return 0;
 737 }