src/tools/eneconv.c

   1 /*
   2  * $Id$
   3  *
   4  *       This source code is part of
   5  *
   6  *        G   R   O   M   A   C   S
   7  *
   8  * GROningen MAchine for Chemical Simulations
   9  *
  10  *               VERSION 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * GRowing Old MAkes el Chrono Sweat
  28  */
  29 static char *SRCID_eneconv_c = "$Id$";
  30
  31 #include <string.h>
  32 #include <math.h>
  33 #include "string2.h"
  34 #include "typedefs.h"
  35 #include "smalloc.h"
  36 #include "statutil.h"
  37 #include "disre.h"
  38 #include "names.h"
  39 #include "copyrite.h"
  40 #include "macros.h"
  41 #include "enxio.h"
  42
  43 #define TIME_EXPLICIT 0
  44 #define TIME_CONTINUE 1
  45 #define TIME_LAST     2
  46 #ifndef FLT_MAX
  47 #define FLT_MAX 1e36
  48 #endif
  49
  50 bool bRmod(double a,double b)
  51 {
  52   int iq;
  53   double tol = 1e-6;
  54
  55   iq = ((1.0+tol)*a)/b;
  56
  57   if (fabs(a-b*iq) <= tol*a)
  58     return TRUE;
  59   else
  60     return FALSE;
  61 }
  62
  63 static bool same_time(real t1,real t2)
  64 {
  65   const real tol=1e-5;
  66
  67   return (fabs(t1-t2) < tol);
  68 }
  69
  70
  71 bool bRgt(double a,double b)
  72 {
  73   double tol = 1e-6;
  74
  75   if ( a > (b - tol*(a+b)) )
  76     return TRUE;
  77   else
  78     return FALSE;
  79 }
  80
  81 static void sort_files(char **fnms,real *settime,int nfile)
  82 {
  83     int i,j,minidx;
  84     real timeswap;
  85     char *chptr;
  86
  87     for(i=0;i<nfile;i++) {
  88         minidx=i;
  89         for(j=i+1;j<nfile;j++) {
  90             if(settime[j]<settime[minidx])
  91                 minidx=j;
  92         }
  93         if(minidx!=i) {
  94             timeswap=settime[i];
  95             settime[i]=settime[minidx];
  96             settime[minidx]=timeswap;
  97             chptr=fnms[i];
  98             fnms[i]=fnms[minidx];
  99             fnms[minidx]=chptr;
 100         }
 101     }
 102 }
 103
 104
 105 static int scan_ene_files(char **fnms,int nfiles,real *readtime, real *timestep)
 106 {
 107     /* Check number of energy terms and start time of all files */
 108     int i,in,nre,ndr,nresav=0,step;
 109     real t1,t2;
 110     char      **enm=NULL;
 111     t_drblock *dr=NULL;
 112     t_energy  *ee=NULL;
 113
 114     for(i=0;i<nfiles;i++) {
 115         in = open_enx(fnms[i],"r");
 116         do_enxnms(in,&nre,&enm);
 117
 118         if (i == 0) {
 119             nresav = nre;
 120             snew(ee,nre);
 121             do_enx(in,&t1,&step,&nre,ee,&ndr,dr);
 122             do_enx(in,&t2,&step,&nre,ee,&ndr,dr);
 123             *timestep=t2-t1;
 124             readtime[i]=t1;
 125             close_enx(in);
 126         }
 127         else if (nre != nresav) {
 128             fatal_error(0,"Energy files don't match, different number"
 129                         " of energies (%s)",fnms[i]);
 130         }
 131         else {
 132             do_enx(in,&t1,&step,&nre,ee,&ndr,dr);
 133             readtime[i]=t1;
 134             close_enx(in);
 135         }
 136         fprintf(stderr,"\n");
 137     }
 138     sfree(ee);
 139     return nre;
 140 }
 141
 142
 143 static void edit_files(char **fnms,int nfiles,real *readtime, real
 144                        *settime, int *cont_type, bool bSetTime,bool bSort)
 145 {
 146     int i;
 147     bool ok;
 148     char inputstring[STRLEN],*chptr;
 149
 150     if(bSetTime) {
 151       if(nfiles==1)
 152         fprintf(stderr,"\n\nEnter the new start time:\n\n");
 153       else
 154         fprintf(stderr,"\n\nEnter the new start time for each file.\n"
 155                 "There are two special options, both disables sorting:\n\n"
 156                 "c (continue) - The start time is taken from the end\n"
 157                 "of the previous file. Use it when your continuation run\n"
 158                 "restarts with t=0 and there is no overlap.\n\n"
 159                 "l (last) - The time in this file will be changed the\n"
 160                 "same amount as in the previous. Use it when the time in the\n"
 161                 "new run continues from the end of the previous one,\n"
 162                 "since this takes possible overlap into account.\n\n");
 163
 164       fprintf(stderr,"          File             Current start       New start\n"
 165               "---------------------------------------------------------\n");
 166
 167       for(i=0;i<nfiles;i++) {
 168         fprintf(stderr,"%25s   %10.3f             ",fnms[i],readtime[i]);
 169         ok=FALSE;
 170         do {
 171           fgets(inputstring,STRLEN-1,stdin);
 172           inputstring[strlen(inputstring)-1]=0;
 173
 174           if(inputstring[0]=='c' || inputstring[0]=='C') {
 175             cont_type[i]=TIME_CONTINUE;
 176             bSort=FALSE;
 177             ok=TRUE;
 178             settime[i]=FLT_MAX;
 179           }
 180           else if(inputstring[0]=='l' ||
 181                   inputstring[0]=='L') {
 182             cont_type[i]=TIME_LAST;
 183             bSort=FALSE;
 184             ok=TRUE;
 185             settime[i]=FLT_MAX;
 186           }
 187           else {
 188             settime[i]=strtod(inputstring,&chptr);
 189             if(chptr==inputstring) {
 190               fprintf(stderr,"Try that again: ");
 191             }
 192             else {
 193               cont_type[i]=TIME_EXPLICIT;
 194               ok=TRUE;
 195             }
 196           }
 197         } while (!ok);
 198       }
 199       if(cont_type[0]!=TIME_EXPLICIT) {
 200         cont_type[0]=TIME_EXPLICIT;
 201         settime[0]=0;
 202       }
 203     }
 204     else
 205         for(i=0;i<nfiles;i++)
 206             settime[i]=readtime[i];
 207
 208     if(bSort && (nfiles>1))
 209         sort_files(fnms,settime,nfiles);
 210     else
 211         fprintf(stderr,"Sorting disabled.\n");
 212
 213
 214     /* Write out the new order and start times */
 215     fprintf(stderr,"\nSummary of files and start times used:\n\n"
 216             "          File                Start time\n"
 217             "-----------------------------------------\n");
 218     for(i=0;i<nfiles;i++)
 219         switch(cont_type[i]) {
 220         case TIME_EXPLICIT:
 221             fprintf(stderr,"%25s   %10.3f\n",fnms[i],settime[i]);
 222             break;
 223         case TIME_CONTINUE:
 224             fprintf(stderr,"%25s        Continue from end of last file\n",fnms[i]);
 225             break;
 226         case TIME_LAST:
 227             fprintf(stderr,"%25s        Change by same amount as last file\n",fnms[i]);
 228             break;
 229         }
 230     fprintf(stderr,"\n");
 231
 232     settime[nfiles]=FLT_MAX;
 233     cont_type[nfiles]=TIME_EXPLICIT;
 234     readtime[nfiles]=FLT_MAX;
 235 }
 236
 237
 238 static void copy_ee(t_energy *src, t_energy *dst, int nre)
 239 {
 240   int i;
 241
 242   for(i=0;i<nre;i++) {
 243     dst[i].e=src[i].e;
 244     dst[i].esum=src[i].esum;
 245     dst[i].eav=src[i].eav;
 246   }
 247 }
 248
 249
 250 static void remove_last_eeframe(t_energy *lastee, int laststep,
 251                                 t_energy *ee, int nre)
 252 {
 253     int i;
 254     int p=laststep+1;
 255     double sigmacorr;
 256
 257     for(i=0;i<nre;i++) {
 258         lastee[i].esum-=ee[i].e;
 259         sigmacorr=lastee[i].esum-(p-1)*ee[i].e;
 260         lastee[i].eav-=(sigmacorr*sigmacorr)/((p-1)*p);
 261     }
 262 }
 263
 264
 265
 266 static void update_ee(t_energy *lastee,int laststep,
 267                       t_energy *startee,int startstep,
 268                       t_energy *ee, int step,
 269                       t_energy *outee, int nre)
 270 {
 271   int i;
 272   double sigmacorr,nom,denom;
 273   double prestart_esum;
 274   double prestart_sigma;
 275
 276   for(i=0;i<nre;i++) {
 277           outee[i].e=ee[i].e;
 278       /* add statistics from earlier file if present */
 279       if(laststep>0) {
 280           outee[i].esum=lastee[i].esum+ee[i].esum;
 281           nom=(lastee[i].esum*(step+1)-ee[i].esum*(laststep));
 282           denom=((step+1.0)*(laststep)*(step+1.0+laststep));
 283           sigmacorr=nom*nom/denom;
 284           outee[i].eav=lastee[i].eav+ee[i].eav+sigmacorr;
 285       }
 286       else {
 287           /* otherwise just copy to output */
 288           outee[i].esum=ee[i].esum;
 289           outee[i].eav=ee[i].eav;
 290       }
 291
 292       /* if we didnt start to write at the first frame
 293        * we must compensate the statistics for this
 294        * there are laststep frames in the earlier file
 295        * and step+1 frames in this one.
 296        */
 297     if(startstep>0) {
 298       int q=laststep+step;
 299       int p=startstep+1;
 300       prestart_esum=startee[i].esum-startee[i].e;
 301       sigmacorr=prestart_esum-(p-1)*startee[i].e;
 302       prestart_sigma=startee[i].eav-
 303         sigmacorr*sigmacorr/(p*(p-1));
 304       sigmacorr=prestart_esum/(p-1)-
 305         outee[i].esum/(q);
 306       outee[i].esum-=prestart_esum;
 307       outee[i].eav=outee[i].eav-prestart_sigma-
 308         sigmacorr*sigmacorr*((p-1)*q)/(q-p+1);
 309     }
 310
 311     if((outee[i].eav/(laststep+step+1))<(GMX_REAL_EPS))
 312       outee[i].eav=0;
 313   }
 314 }
 315
 316
 317 static void update_last_ee(t_energy *lastee, int laststep,
 318                            t_energy *ee,int step,int nre)
 319 {
 320     t_energy *tmp;
 321     snew(tmp,nre);
 322     update_ee(lastee,laststep,NULL,0,ee,step,tmp,nre);
 323     copy_ee(tmp,lastee,nre);
 324     sfree(tmp);
 325 }
 326
 327
 328 int main(int argc,char *argv[])
 329 {
 330   static char *desc[] = {
 331     "When [TT]-f[tt] is [IT]not[it] specified:[BR]",
 332     "Concatenates several energy files in sorted order.",
 333     "In case of double time frames the one",
 334     "in the later file is used. By specifying [TT]-settime[tt] you will be",
 335     "asked for the start time of each file. The input files are taken",
 336     "from the command line,",
 337     "such that the command [TT]eneconv -o fixed.edr *.edr[tt] should do",
 338     "the trick. [PAR]",
 339     "With [TT]-f[tt] specified:[BR]",
 340     "Reads one energy file and writes another, applying the [TT]-dt[tt],",
 341     "[TT]-offset[tt], [TT]-t0[tt] and [TT]-settime[tt] options and",
 342     "converting to a different format if necessary (indicated by file",
 343     "extentions).[PAR]",
 344     "[TT]-settime[tt] is applied first, then [TT]-dt[tt]/[TT]-offset[tt]",
 345     "followed by [TT]-b[tt] and [TT]-e[tt] to select which frames to write."
 346   };
 347
 348   int       in,out=0;
 349   t_energy  *ee,*lastee,*outee,*startee;
 350   int       step,laststep,outstep,startstep;
 351   int       nre,nfile,i,j,ndr;
 352   real      t=0,outt=-1;
 353   char      **fnms;
 354   char      **enm=NULL;
 355   t_drblock *dr=NULL;
 356   real      *readtime,*settime,timestep,t1,tadjust;
 357   char      inputstring[STRLEN],*chptr;
 358   bool      ok;
 359   int       *cont_type;
 360   bool      bNewFile,bFirst,bNewOutput;
 361
 362   t_filenm fnm[] = {
 363     { efENX, "-f", NULL,    ffREAD  },
 364     { efENX, "-o", "fixed", ffOPTWR },
 365   };
 366
 367 #define NFILE asize(fnm)
 368
 369   static real  delta_t=0.0, toffset=0;
 370   static bool  bSetTime=FALSE;
 371   static bool  bSort=TRUE;
 372   static real  begin=-1;
 373   static real  end=-1;
 374
 375   t_pargs pa[] = {
 376       { "-b",        FALSE, etREAL, {&begin},
 377         "First time to use"},
 378       { "-e",        FALSE, etREAL, {&end},
 379         "Last time to use"},
 380       { "-dt",       FALSE, etREAL, {&delta_t},
 381         "Only write out frame when t MOD dt = offset" },
 382       { "-offset",   FALSE, etREAL, {&toffset},
 383         "Time offset for -dt option" },
 384       { "-settime",  FALSE, etBOOL, {&bSetTime},
 385         "Change starting time interactively" },
 386       { "-sort",     FALSE, etBOOL, {&bSort},
 387         "Sort energy files (not frames)"}
 388   };
 389
 390   CopyRight(stderr,argv[0]);
 391   parse_common_args(&argc,argv,PCA_NOEXIT_ON_ARGS,TRUE,
 392                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL);
 393   tadjust=0;
 394   snew(fnms,argc);
 395   nfile=0;
 396   outstep=laststep=startstep=0;
 397
 398   for(i=1; (i<argc); i++)
 399         fnms[nfile++]=argv[i];
 400   if(nfile==0)
 401     nfile=1;
 402   snew(settime,nfile+1);
 403   snew(readtime,nfile+1);
 404   snew(cont_type,nfile+1);
 405
 406   if (!opt2bSet("-f",NFILE,fnm)) {
 407     if(!nfile)
 408           fatal_error(0,"No input files!");
 409  }
 410   else {
 411     /* get the single filename */
 412     fnms[0]=opt2fn("-f",NFILE,fnm);
 413   }
 414
 415   nre=scan_ene_files(fnms,nfile,readtime,&timestep);
 416   edit_files(fnms,nfile,readtime,settime,cont_type,bSetTime,bSort);
 417
 418   snew(ee,nre);
 419   snew(outee,nre);
 420
 421   if(nfile>1)
 422     snew(lastee,nre);
 423   else
 424     lastee=NULL;
 425
 426   if(begin>0)
 427     snew(startee,nre);
 428   else
 429     startee=NULL;
 430
 431   bFirst=TRUE;
 432
 433   for(i=0;i<nfile;i++) {
 434       bNewFile=TRUE;
 435       bNewOutput=TRUE;
 436     in=open_enx(fnms[i],"r");
 437     do_enxnms(in,&nre,&enm);
 438     if(i==0) {
 439         /* write names to the output file */
 440         out=open_enx(opt2fn("-o",NFILE,fnm),"w");
 441         do_enxnms(out,&nre,&enm);
 442     }
 443
 444     /* start reading from the next file */
 445     while((t<(settime[i+1]-GMX_REAL_EPS)) &&
 446           do_enx(in,&t1,&step,&nre,ee,&ndr,dr)) {
 447       if(bNewFile) {
 448           tadjust=settime[i]-t1;
 449           if(cont_type[i+1]==TIME_LAST) {
 450               settime[i+1]=readtime[i+1]-readtime[i]+settime[i];
 451               cont_type[i+1]=TIME_EXPLICIT;
 452           }
 453         bNewFile=FALSE;
 454       }
 455       t=tadjust+t1;
 456
 457       if((end>0) && (t>(end+GMX_REAL_EPS))) {
 458           i=nfile;
 459           break;
 460       }
 461
 462       if(t>=(begin-GMX_REAL_EPS)) {
 463           if((bFirst)) {
 464               bFirst=FALSE;
 465               if(startee!=NULL)
 466                   copy_ee(ee,startee,nre);
 467               startstep=step;
 468           }
 469           update_ee(lastee,laststep,startee,startstep,ee,step,outee,nre);
 470           outstep=laststep+step-startstep;
 471       }
 472
 473       /* determine if we should write it */
 474       if((t<(settime[i+1]-GMX_REAL_EPS)) && (t>=begin) &&
 475          ((delta_t==0) || (bRmod(t-toffset,delta_t)))) {
 476           outt=t;
 477           if(bNewOutput) {
 478               bNewOutput=FALSE;
 479               fprintf(stderr,"\nContinue writing frames from t=%g, step=%d\n",
 480                       t,outstep);
 481           }
 482           do_enx(out,&outt,&outstep,&nre,outee,&ndr,dr);
 483           fprintf(stderr,"\rWriting step %d, time %f        ",outstep,outt);
 484       }
 485     }
 486     /* copy statistics to old */
 487     if(lastee!=NULL) {
 488         update_last_ee(lastee,laststep,ee,step,nre);
 489         laststep+=step;
 490         /* remove the last frame from statistics since gromacs2.0 repeats it in the next file */
 491         remove_last_eeframe(lastee,laststep,ee,nre);
 492         /* the old part now has (laststep) values, and the new (step+1) */
 493         printf("laststep=%d step=%d\n",laststep,step);
 494     }
 495
 496     /* set the next time from the last in previous file */
 497     if(cont_type[i+1]==TIME_CONTINUE) {
 498         settime[i+1]=outt;
 499         /* in this case we have already written the last frame of
 500          * previous file, so update begin to avoid doubling it
 501          * with the start of the next file
 502          */
 503         begin=outt+0.5*timestep;
 504         /* cont_type[i+1]==TIME_EXPLICIT; */
 505     }
 506
 507     if((outt<end) && (i<(nfile-1)) &&
 508        (outt<(settime[i+1]-1.5*timestep)))
 509       fprintf(stderr,
 510               "\nWARNING: There might be a gap around t=%g\n",t);
 511
 512     /* move energies to lastee */
 513     close_enx(in);
 514
 515     fprintf(stderr,"\n");
 516   }
 517   if(outstep==0)
 518       fprintf(stderr,"No frames written.\n");
 519   else
 520       fprintf(stderr,"Last frame written was at step %d, time %f\n",outstep,outt);
 521
 522   thanx(stdout);
 523   return 0;
 524 }