src/tools/trjcat.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  * Great Red Oystrich Makes All Chemists Sane
  28  */
  29 static char *SRCID_trjcat_c = "$Id$";
  30
  31 #include <string.h>
  32 #include <math.h>
  33 #include <unistd.h>
  34 #include "macros.h"
  35 #include "sysstuff.h"
  36 #include "smalloc.h"
  37 #include "typedefs.h"
  38 #include "copyrite.h"
  39 #include "gmxfio.h"
  40 #include "tpxio.h"
  41 #include "trnio.h"
  42 #include "statutil.h"
  43 #include "futil.h"
  44 #include "pdbio.h"
  45 #include "confio.h"
  46 #include "names.h"
  47 #include "rdgroup.h"
  48 #include "vec.h"
  49 #include "xtcio.h"
  50 #include "do_fit.h"
  51 #include "rmpbc.h"
  52 #include "wgms.h"
  53 #include "magic.h"
  54 #include "binio.h"
  55 #include "pbc.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 void scan_trj_files(char **fnms,int nfiles,real *readtime, real *timestep)
  65 {
  66     /* Check start time of all files */
  67     int i,status,oldnatoms=0,natoms,step;
  68     real t,t2;
  69     rvec *x;
  70     matrix box;
  71
  72
  73     for(i=0;i<nfiles;i++) {
  74       natoms=read_first_x(&status,fnms[i],&t,&x,box);
  75       if(i==0) {
  76         oldnatoms=natoms;
  77         read_next_x(status,&t2,natoms,x,  box);
  78         *timestep=t2-t;
  79         sfree(x);
  80       }
  81       else {
  82         if(natoms!=oldnatoms)
  83           fatal_error(0,"Different number of atoms in files");
  84       }
  85       readtime[i]=t;
  86       close_trj(status);
  87     }
  88     fprintf(stderr,"\n");
  89 }
  90
  91
  92 static void sort_files(char **fnms,real *settime,int nfile)
  93 {
  94     int i,j,minidx;
  95     real timeswap;
  96     char *chptr;
  97
  98     for(i=0;i<nfile;i++) {
  99         minidx=i;
 100         for(j=i+1;j<nfile;j++) {
 101             if(settime[j]<settime[minidx])
 102                 minidx=j;
 103         }
 104         if(minidx!=i) {
 105             timeswap=settime[i];
 106             settime[i]=settime[minidx];
 107             settime[minidx]=timeswap;
 108             chptr=fnms[i];
 109             fnms[i]=fnms[minidx];
 110             fnms[minidx]=chptr;
 111         }
 112     }
 113 }
 114
 115
 116 static void edit_files(char **fnms,int nfiles,real *readtime, real
 117                        *settime, int *cont_type, bool bSetTime,bool bSort)
 118 {
 119     int i;
 120     bool ok;
 121     char inputstring[STRLEN],*chptr;
 122
 123     if(bSetTime) {
 124         fprintf(stderr,"\n\nEnter the new start time for each file.\n"
 125                 "There are two special options, both disable sorting:\n\n"
 126                 "c (continue) - The start time is taken from the end\n"
 127                 "of the previous file. Use it when your continuation run\n"
 128                 "restarts with t=0.\n\n"
 129                 "l (last) - The time in this file will be changed the\n"
 130                 "same amount as in the previous. Use it when the time in the\n"
 131                 "new run continues from the end of the previous one,\n"
 132                 "since this takes possible overlap into account.\n\n");
 133
 134           fprintf(stderr,"          File             Current start       New start\n"
 135                   "---------------------------------------------------------\n");
 136
 137           for(i=0;i<nfiles;i++) {
 138               fprintf(stderr,"%25s   %10.3f             ",fnms[i],readtime[i]);
 139               ok=FALSE;
 140               do {
 141                   fgets(inputstring,STRLEN-1,stdin);
 142                   inputstring[strlen(inputstring)-1]=0;
 143
 144                   if(inputstring[0]=='c' || inputstring[0]=='C') {
 145                     cont_type[i]=TIME_CONTINUE;
 146                     bSort=FALSE;
 147                     ok=TRUE;
 148                     settime[i]=FLT_MAX;
 149                   }
 150                   else if(inputstring[0]=='l' ||
 151                           inputstring[0]=='L') {
 152                     cont_type[i]=TIME_LAST;
 153                     bSort=FALSE;
 154                     ok=TRUE;
 155                     settime[i]=FLT_MAX;
 156                   }
 157                   else {
 158                     settime[i]=strtod(inputstring,&chptr);
 159                     if(chptr==inputstring) {
 160                       fprintf(stderr,"'%s' not recognized as a floating point number, 'c' or 'l'. "
 161                               "Try again: ",inputstring);
 162                     }
 163                     else {
 164                       cont_type[i]=TIME_EXPLICIT;
 165                       ok=TRUE;
 166                     }
 167                   }
 168               } while (!ok);
 169           }
 170           if(cont_type[0]!=TIME_EXPLICIT) {
 171               cont_type[0]=TIME_EXPLICIT;
 172               settime[0]=0;
 173           }
 174     }
 175     else
 176         for(i=0;i<nfiles;i++)
 177             settime[i]=readtime[i];
 178
 179     if(!bSort)
 180         fprintf(stderr,"Sorting disabled.\n");
 181     else
 182         sort_files(fnms,settime,nfiles);
 183
 184
 185     /* Write out the new order and start times */
 186     fprintf(stderr,"\nSummary of files and start times used:\n\n"
 187             "          File                Start time\n"
 188             "-----------------------------------------\n");
 189     for(i=0;i<nfiles;i++)
 190         switch(cont_type[i]) {
 191         case TIME_EXPLICIT:
 192             fprintf(stderr,"%25s   %10.3f\n",fnms[i],settime[i]);
 193             break;
 194         case TIME_CONTINUE:
 195             fprintf(stderr,"%25s        Continue from last file\n",fnms[i]);
 196             break;
 197         case TIME_LAST:
 198             fprintf(stderr,"%25s        Change by same amount as last file\n",fnms[i]);
 199             break;
 200         }
 201     fprintf(stderr,"\n");
 202
 203     settime[nfiles]=FLT_MAX;
 204     cont_type[nfiles]=TIME_EXPLICIT;
 205     readtime[nfiles]=FLT_MAX;
 206 }
 207
 208
 209 bool bRmod(double a,double b)
 210 {
 211   int iq;
 212   double tol = 1e-6;
 213
 214   iq = ((1.0+tol)*a)/b;
 215
 216   if (fabs(a-b*iq) <= tol*a)
 217     return TRUE;
 218   else
 219     return FALSE;
 220 }
 221
 222 int main(int argc,char *argv[])
 223 {
 224   static char *desc[] = {
 225       "trjcat concatenates several input trajectory files in sorted order. ",
 226       "In case of double time frames the one in the later file is used. ",
 227       "By specifying [TT]-settime[tt] you will be asked for the start time ",
 228       "of each file. The input files are taken from the command line, ",
 229       "such that a command like [TT]trjconv -o fixed.trr *.trr[tt] should do ",
 230       "the trick."
 231   };
 232   static bool  bVels=TRUE;
 233   static int   prec=3;
 234   static bool  bSort=TRUE;
 235   static bool  bSetTime=FALSE;
 236   static real  begin=-1;
 237   static real  end=-1;
 238
 239   t_pargs pa[] = {
 240       { "-b",        FALSE, etREAL, {&begin},
 241         "First time to use"},
 242       { "-e",        FALSE, etREAL, {&end},
 243         "Last time to use"},
 244     { "-prec", FALSE,  etINT,  {&prec},
 245       "Precision for .xtc and .gro writing in number of decimal places" },
 246     { "-vel", FALSE, etBOOL, {&bVels},
 247       "Read and write velocities if possible" },
 248     { "-settime",  FALSE, etBOOL, {&bSetTime},
 249       "Change starting time interactively" },
 250     { "-sort",     FALSE, etBOOL, {&bSort},
 251       "Sort trajectory files (not frames)"}
 252   };
 253
 254   FILE         *out=NULL;
 255   int          status,ftp,ftpin,i,frame,natoms=0,step,trjout=0;
 256   rvec         *x,*v;
 257   real         xtcpr,t,t0=-1,t1;
 258   matrix       box;
 259   t_topology   top;
 260   char         **fnms;
 261   bool         bNewFile,bHaveX=FALSE,bHaveV=FALSE;
 262   char         *in_file,*out_file;
 263   int          xdr=0,earliersteps,nfile,*cont_type,last_ok_step;
 264   real         *readtime,*settime,tstart,last_ok_t=-1,timestep;
 265
 266   t_filenm fnm[] = {
 267       { efTRX, "-o", "trajout", ffWRITE }
 268   };
 269
 270 #define NFILE asize(fnm)
 271
 272   CopyRight(stderr,argv[0]);
 273   parse_common_args(&argc,argv,PCA_NOEXIT_ON_ARGS,TRUE,
 274                     NFILE,fnm,asize(pa),pa,asize(desc),desc,
 275                     0,NULL);
 276
 277   /* prec is in nr of decimal places, xtcprec is a multiplication factor: */
 278   xtcpr=1;
 279   for (i=0; i<prec; i++)
 280       xtcpr*=10;
 281
 282   nfile=0;
 283
 284   snew(fnms,argc);
 285   for(i=1; (i<argc); i++)
 286       fnms[nfile++]=argv[i];
 287
 288   if(!nfile)
 289       fatal_error(0,"No input files!");
 290
 291   snew(settime,nfile+1);
 292   snew(readtime,nfile+1);
 293   snew(cont_type,nfile+1);
 294
 295   scan_trj_files(fnms,nfile,readtime,&timestep);
 296   edit_files(fnms,nfile,readtime,settime,cont_type,bSetTime,bSort);
 297
 298   t=-1;
 299   earliersteps=0;
 300   out_file=opt2fn("-o",NFILE,fnm);
 301   ftp=fn2ftp(out_file);
 302
 303   bVels= ((ftp==efTRR) ||(ftp==efTRJ) || (ftp==efGRO));
 304   /* Not checking input format, could be dangerous :-) */
 305   if (!bVels) {
 306     bHaveX=TRUE;
 307     bHaveV=FALSE;
 308     v=NULL;
 309   }
 310
 311
 312   frame=-1;
 313   /* Lets stitch up some files */
 314   for(i=0;i<nfile;i++) {
 315       /* Open next file */
 316       if(bVels)
 317           natoms=read_first_x_or_v(&status,fnms[i],&tstart,&x,&v,box);
 318       else
 319           natoms=read_first_x(&status,fnms[i],&tstart,&x,box);
 320       if(cont_type[i]==TIME_EXPLICIT)
 321           t0=settime[i]-tstart;
 322       t1=tstart;
 323
 324       bNewFile=TRUE;
 325       if(i==0) {
 326           switch(ftp) {
 327           case efXTC:
 328               xdr = open_xtc(out_file,"w");
 329               break;
 330           case efG87:
 331               out=ffopen(out_file,"w");
 332               break;
 333           case efTRR:
 334           case efTRJ:
 335               out=NULL;
 336               trjout = open_tpx(out_file,"w");
 337               break;
 338           case efGRO:
 339           case efG96:
 340           case efPDB:
 341                   out=ffopen(out_file,"w");
 342               break;
 343           }
 344       }
 345
 346       do {
 347           /* set the new time */
 348           t=t0+t1;
 349           if((end>0) && (t>(end+GMX_REAL_EPS))) {
 350               i=nfile;
 351               break;
 352           }
 353           /* determine if we should write it */
 354           if((t<(settime[i+1]-0.5*timestep)) && (t>=begin)) {
 355               frame++;
 356               last_ok_t=t;
 357               if(bNewFile) {
 358                   fprintf(stderr,"\nContinue writing frames from t=%g, frame=%d      \n",t,frame);
 359                   bNewFile=FALSE;
 360               }
 361
 362               switch(ftp) {
 363               case efTRJ:
 364               case efTRR:
 365                   fwrite_trn(trjout,frame,t,0,box,
 366                              natoms,x,v,NULL);
 367                   break;
 368               case efXTC:
 369                   write_xtc(xdr,natoms,frame,t,box,x,xtcpr);
 370                   break;
 371               default:
 372                   fatal_error(0,"DHE, ftp=%d\n",ftp);
 373               }
 374
 375               if ( ((frame % 10) == 0) || (frame < 10) )
 376                 fprintf(stderr," ->  frame %6d time %8.3f      \r",frame,t);
 377           }
 378       } while((t<settime[i+1]) &&
 379               ((bVels && read_next_x_or_v(status,&t1,natoms,x,v,box)) ||
 380                (read_next_x(status,&t1,natoms,x,  box))));
 381
 382       close_trj(status);
 383
 384       earliersteps+=step;
 385
 386       /* set the next time from the last in previous file */
 387       if(cont_type[i+1]==TIME_CONTINUE) {
 388           begin=t+0.5*timestep;
 389           settime[i+1]=t;
 390           cont_type[i+1]=TIME_EXPLICIT;
 391       }
 392       else if(cont_type[i+1]==TIME_LAST)
 393           begin=t+0.5*timestep;
 394
 395       if(cont_type[i+1]==TIME_EXPLICIT)
 396         if((i<(nfile-1)) &&
 397            (t<(settime[i+1]-1.5*timestep)))
 398           fprintf(stderr,
 399                   "WARNING: Frames around t=%f have a different spacing than the rest,\n"
 400                   "might be a gap or overlap that couldn't be corrected automatically.\n",t);
 401   }
 402
 403   if (ftp == efXTC)
 404       close_xtc(xdr);
 405   else if (out != NULL)
 406       fclose(out);
 407   fprintf(stderr,"\nLast frame written was %d, time %f\n",frame,last_ok_t);
 408
 409   thanx(stdout);
 410
 411   return 0;
 412 }