src/gmxlib/trxio.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 forxd
  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  * GROningen Mixture of Alchemy and Childrens' Stories
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <ctype.h>
  40 #include "sysstuff.h"
  41 #include "string2.h"
  42 #include "smalloc.h"
  43 #include "pbc.h"
  44 #include "statutil.h"
  45 #include "gmxfio.h"
  46 #include "trnio.h"
  47 #include "names.h"
  48 #include "vec.h"
  49 #include "futil.h"
  50 #include "gmxfio.h"
  51 #include "xtcio.h"
  52 #include "pdbio.h"
  53 #include "confio.h"
  54 #include "checkpoint.h"
  55 #include "wgms.h"
  56 #include <math.h>
  57
  58 /* defines for frame counter output */
  59 static int __frame=NOTSET;
  60 #define SKIP1   10
  61 #define SKIP2  100
  62 #define SKIP3 1000
  63 #define INITCOUNT __frame=-1
  64
  65
  66 /* frames for read_first/next_x */
  67 static t_trxframe *xframe=NULL;
  68 static int nxframe=0;
  69
  70
  71 int nframes_read(void)
  72 {
  73   return __frame;
  74 }
  75
  76 static void printcount_(const output_env_t oenv,const char *l,real t)
  77 {
  78   if ((__frame < 2*SKIP1 || __frame % SKIP1 == 0) &&
  79       (__frame < 2*SKIP2 || __frame % SKIP2 == 0) &&
  80       (__frame < 2*SKIP3 || __frame % SKIP3 == 0))
  81     fprintf(stderr,"\r%-14s %6d time %8.3f   ",l,__frame,conv_time(oenv,t));
  82 }
  83
  84 static void printcount(const output_env_t oenv,real t,bool bSkip)
  85 {
  86   __frame++;
  87   printcount_(oenv,bSkip ? "Skipping frame" : "Reading frame",t);
  88 }
  89
  90 static void printlast(const output_env_t oenv,real t)
  91 {
  92   printcount_(oenv,"Last frame",t);
  93   fprintf(stderr,"\n");
  94 }
  95
  96 static void printincomp(t_trxframe *fr)
  97 {
  98   if (fr->not_ok & HEADER_NOT_OK)
  99     fprintf(stderr,"WARNING: Incomplete header: nr %d time %g\n",
 100             __frame+1,fr->time);
 101   else if (fr->not_ok)
 102     fprintf(stderr,"WARNING: Incomplete frame: nr %d time %g\n",
 103             __frame+1,fr->time);
 104 }
 105
 106 int prec2ndec(real prec)
 107 {
 108   if (prec <= 0)
 109     gmx_fatal(FARGS,"DEATH HORROR prec (%g) <= 0 in prec2ndec",prec);
 110
 111   return (int)(log(prec)/log(10.0)+0.5);
 112 }
 113
 114 /* Globals for gromos-87 input */
 115 typedef enum { effXYZ, effXYZBox, effG87, effG87Box, effNR } eFileFormat;
 116 static eFileFormat eFF;
 117 static int         NATOMS;
 118 static double      DT,BOX[3];
 119 static bool        bReadBox;
 120
 121 void clear_trxframe(t_trxframe *fr,bool bFirst)
 122 {
 123   fr->not_ok  = 0;
 124   fr->bTitle  = FALSE;
 125   fr->bStep   = FALSE;
 126   fr->bTime   = FALSE;
 127   fr->bLambda = FALSE;
 128   fr->bAtoms  = FALSE;
 129   fr->bPrec   = FALSE;
 130   fr->bX      = FALSE;
 131   fr->bV      = FALSE;
 132   fr->bF      = FALSE;
 133   fr->bBox    = FALSE;
 134   if (bFirst) {
 135     fr->flags  = 0;
 136     fr->bDouble= FALSE;
 137     fr->natoms = -1;
 138     fr->t0     = 0;
 139     fr->tpf    = 0;
 140     fr->tppf   = 0;
 141     fr->title  = NULL;
 142     fr->step   = 0;
 143     fr->time   = 0;
 144     fr->lambda = 0;
 145     fr->atoms  = NULL;
 146     fr->prec   = 0;
 147     fr->x      = NULL;
 148     fr->v      = NULL;
 149     fr->f      = NULL;
 150     clear_mat(fr->box);
 151     fr->bPBC   = FALSE;
 152     fr->ePBC   = -1;
 153   }
 154 }
 155
 156 void set_trxframe_ePBC(t_trxframe *fr,int ePBC)
 157 {
 158   fr->bPBC = (ePBC == -1);
 159   fr->ePBC = ePBC;
 160 }
 161
 162 int write_trxframe_indexed(int fnum,t_trxframe *fr,int nind,atom_id *ind,
 163                            gmx_conect gc)
 164 {
 165   char title[STRLEN];
 166   rvec *xout=NULL,*vout=NULL,*fout=NULL;
 167   int  i;
 168   real prec;
 169
 170   if (fr->bPrec)
 171     prec = fr->prec;
 172   else
 173     prec = 1000.0;
 174
 175   switch (gmx_fio_getftp(fnum)) {
 176   case efTRJ:
 177   case efTRR:
 178     break;
 179   default:
 180     if (!fr->bX)
 181       gmx_fatal(FARGS,"Need coordinates to write a %s trajectory",
 182                   ftp2ext(gmx_fio_getftp(fnum)));
 183     break;
 184   }
 185
 186   switch (gmx_fio_getftp(fnum)) {
 187   case efTRJ:
 188   case efTRR:
 189     if (fr->bV) {
 190       snew(vout,nind);
 191       for(i=0; i<nind; i++)
 192         copy_rvec(fr->v[ind[i]],vout[i]);
 193     }
 194     if (fr->bF) {
 195       snew(fout,nind);
 196       for(i=0; i<nind; i++)
 197         copy_rvec(fr->f[ind[i]],fout[i]);
 198     }
 199   case efXTC:
 200   case efG87:
 201     if (fr->bX) {
 202       snew(xout,nind);
 203       for(i=0; i<nind; i++)
 204         copy_rvec(fr->x[ind[i]],xout[i]);
 205     }
 206     break;
 207   default:
 208     break;
 209   }
 210
 211   switch (gmx_fio_getftp(fnum)) {
 212   case efXTC:
 213     write_xtc(fnum,nind,fr->step,fr->time,fr->box,xout,prec);
 214     break;
 215   case efTRJ:
 216   case efTRR:
 217     fwrite_trn(fnum,nframes_read(),
 218                fr->time,fr->step,fr->box,nind,xout,vout,fout);
 219     break;
 220   case efGRO:
 221   case efPDB:
 222   case efBRK:
 223   case efENT:
 224     if (!fr->bAtoms)
 225       gmx_fatal(FARGS,"Can not write a %s file without atom names",
 226                   ftp2ext(gmx_fio_getftp(fnum)));
 227     sprintf(title,"frame t= %.3f",fr->time);
 228     if (gmx_fio_getftp(fnum) == efGRO)
 229       write_hconf_indexed_p(gmx_fio_getfp(fnum),title,fr->atoms,nind,ind,
 230                             prec2ndec(prec),
 231                             fr->x,fr->bV ? fr->v : NULL,fr->box);
 232     else
 233       write_pdbfile_indexed(gmx_fio_getfp(fnum),title,fr->atoms,
 234                             fr->x,-1,fr->box,0,fr->step,nind,ind,gc);
 235     break;
 236   case efG87:
 237     write_gms(gmx_fio_getfp(fnum),nind,xout,fr->box);
 238     break;
 239   case efG96:
 240     write_g96_conf(gmx_fio_getfp(fnum),fr,nind,ind);
 241     break;
 242   default:
 243     gmx_fatal(FARGS,"Sorry, write_trxframe_indexed can not write %s",
 244                 ftp2ext(gmx_fio_getftp(fnum)));
 245     break;
 246   }
 247
 248   switch (gmx_fio_getftp(fnum)) {
 249   case efTRN:
 250   case efTRJ:
 251   case efTRR:
 252     if (vout) sfree(vout);
 253     if (fout) sfree(fout);
 254   case efXTC:
 255   case efG87:
 256     sfree(xout);
 257     break;
 258   default:
 259     break;
 260   }
 261
 262   return 0;
 263 }
 264
 265 int write_trxframe(int fnum,t_trxframe *fr,gmx_conect gc)
 266 {
 267   char title[STRLEN];
 268   real prec;
 269
 270   if (fr->bPrec)
 271     prec = fr->prec;
 272   else
 273     prec = 1000.0;
 274
 275   switch (gmx_fio_getftp(fnum)) {
 276   case efTRJ:
 277   case efTRR:
 278     break;
 279   default:
 280     if (!fr->bX)
 281       gmx_fatal(FARGS,"Need coordinates to write a %s trajectory",
 282                   ftp2ext(gmx_fio_getftp(fnum)));
 283     break;
 284   }
 285
 286   switch (gmx_fio_getftp(fnum)) {
 287   case efXTC:
 288     write_xtc(fnum,fr->natoms,fr->step,fr->time,fr->box,fr->x,prec);
 289     break;
 290   case efTRJ:
 291   case efTRR:
 292     fwrite_trn(fnum,fr->step,fr->time,fr->lambda,fr->box,fr->natoms,
 293                fr->bX ? fr->x:NULL,fr->bV ? fr->v:NULL ,fr->bF ? fr->f:NULL);
 294     break;
 295   case efGRO:
 296   case efPDB:
 297   case efBRK:
 298   case efENT:
 299     if (!fr->bAtoms)
 300       gmx_fatal(FARGS,"Can not write a %s file without atom names",
 301                   ftp2ext(gmx_fio_getftp(fnum)));
 302     sprintf(title,"frame t= %.3f",fr->time);
 303     if (gmx_fio_getftp(fnum) == efGRO)
 304       write_hconf_p(gmx_fio_getfp(fnum),title,fr->atoms,
 305                     prec2ndec(prec),fr->x,fr->bV ? fr->v : NULL,fr->box);
 306     else
 307       write_pdbfile(gmx_fio_getfp(fnum),title,
 308                     fr->atoms,fr->x,fr->bPBC ? fr->ePBC : -1,fr->box,
 309                     0,fr->step,gc);
 310     break;
 311   case efG87:
 312     write_gms(gmx_fio_getfp(fnum),fr->natoms,fr->x,fr->box);
 313     break;
 314   case efG96:
 315     write_g96_conf(gmx_fio_getfp(fnum),fr,-1,NULL);
 316     break;
 317   default:
 318     gmx_fatal(FARGS,"Sorry, write_trxframe can not write %s",
 319                 ftp2ext(gmx_fio_getftp(fnum)));
 320     break;
 321   }
 322
 323   return 0;
 324 }
 325
 326 int write_trx(int fnum,int nind,atom_id *ind,t_atoms *atoms,
 327               int step,real time,matrix box,rvec x[],rvec *v,
 328               gmx_conect gc)
 329 {
 330   t_trxframe fr;
 331
 332   clear_trxframe(&fr,TRUE);
 333   fr.bStep = TRUE;
 334   fr.step = step;
 335   fr.bTime = TRUE;
 336   fr.time = time;
 337   fr.bAtoms = atoms!=NULL;
 338   fr.atoms = atoms;
 339   fr.bX = TRUE;
 340   fr.x = x;
 341   fr.bV = v!=NULL;
 342   fr.v = v;
 343   fr.bBox = TRUE;
 344   copy_mat(box,fr.box);
 345
 346   return write_trxframe_indexed(fnum,&fr,nind,ind,gc);
 347 }
 348
 349 void close_trx(int status)
 350 {
 351   gmx_fio_close(status);
 352 }
 353
 354 int open_trx(const char *outfile,const char *filemode)
 355 {
 356   if (filemode[0]!='w' && filemode[0]!='a')
 357     gmx_fatal(FARGS,"Sorry, write_trx can only write");
 358
 359   return gmx_fio_open(outfile,filemode);
 360 }
 361
 362 static bool gmx_next_frame(int status,t_trxframe *fr)
 363 {
 364   t_trnheader sh;
 365   bool bOK,bRet;
 366
 367   bRet = FALSE;
 368
 369   if (fread_trnheader(status,&sh,&bOK)) {
 370     fr->bDouble=sh.bDouble;
 371     fr->natoms=sh.natoms;
 372     fr->bStep=TRUE;
 373     fr->step=sh.step;
 374     fr->bTime=TRUE;
 375     fr->time=sh.t;
 376     fr->bLambda = TRUE;
 377     fr->lambda = sh.lambda;
 378     fr->bBox = sh.box_size>0;
 379     if (fr->flags & (TRX_READ_X | TRX_NEED_X)) {
 380       if (fr->x==NULL)
 381         snew(fr->x,sh.natoms);
 382       fr->bX = sh.x_size>0;
 383     }
 384     if (fr->flags & (TRX_READ_V | TRX_NEED_V)) {
 385       if (fr->v==NULL)
 386         snew(fr->v,sh.natoms);
 387       fr->bV = sh.v_size>0;
 388     }
 389     if (fr->flags & (TRX_READ_F | TRX_NEED_F)) {
 390       if (fr->f==NULL)
 391         snew(fr->f,sh.natoms);
 392       fr->bF = sh.f_size>0;
 393     }
 394     if (fread_htrn(status,&sh,fr->box,fr->x,fr->v,fr->f))
 395       bRet = TRUE;
 396     else
 397       fr->not_ok = DATA_NOT_OK;
 398   } else
 399     if (!bOK)
 400       fr->not_ok = HEADER_NOT_OK;
 401
 402   return bRet;
 403 }
 404
 405 static void choose_ff(FILE *status)
 406 {
 407   int i,m,c;
 408   int rc;
 409
 410   printf("\n\n");
 411   printf("   Select File Format\n");
 412   printf("---------------------------\n");
 413   printf("1. XYZ File\n");
 414   printf("2. XYZ File with Box\n");
 415   printf("3. Gromos-87 Ascii Trajectory\n");
 416   printf("4. Gromos-87 Ascii Trajectory with Box\n");
 417
 418   do {
 419     printf("\nChoice: ");
 420     fflush(stdout);
 421     do
 422     {
 423       rc = scanf("%d",&i);
 424     }
 425     while (rc!=1);
 426     i--;
 427   } while ((i < 0) || (i >= effNR));
 428   printf("\n");
 429
 430   eFF = (eFileFormat) i;
 431
 432   for(m=0; (m<DIM); m++) BOX[m]=0;
 433
 434   bReadBox = (eFF == effG87Box) || (eFF == effXYZBox);
 435
 436   switch (eFF) {
 437   case effXYZ:
 438   case effXYZBox:
 439     if( 5 != fscanf(status,"%d%lf%lf%lf%lf",&NATOMS,&BOX[XX],&BOX[YY],&BOX[ZZ],&DT))
 440     {
 441       gmx_fatal(FARGS,"Error reading natoms/box in file");
 442     }
 443     break;
 444   case effG87:
 445   case effG87Box:
 446     printf("GROMOS! OH DEAR...\n\n");
 447     printf("Number of atoms ? ");
 448     fflush(stdout);
 449     if (1 != scanf("%d",&NATOMS))
 450     {
 451         gmx_fatal(FARGS,"Error reading natoms in file");
 452     }
 453
 454     printf("Time between timeframes ? ");
 455     fflush(stdout);
 456     if( 1 != scanf("%lf",&DT))
 457     {
 458         gmx_fatal(FARGS,"Error reading dt from file");
 459     }
 460
 461     if (eFF == effG87) {
 462       printf("Box X Y Z ? ");
 463       fflush(stdout);
 464       if(3 != scanf("%lf%lf%lf",&BOX[XX],&BOX[YY],&BOX[ZZ]))
 465       {
 466           gmx_fatal(FARGS,"Error reading box in file");
 467       }
 468     }
 469     do {
 470       c=fgetc(status);
 471       printf("%c",c);
 472     } while (c != '\n');
 473     printf("\n");
 474     fflush(stdout);
 475     break;
 476   default:
 477     printf("Hellow World\n");
 478   }
 479 }
 480
 481 static bool do_read_xyz(FILE *status,int natoms,rvec x[],matrix box)
 482 {
 483   int    i,m;
 484   double x0;
 485
 486   for(i=0; (i<natoms); i++) {
 487     for(m=0; (m<DIM); m++) {
 488       if (fscanf(status,"%lf",&x0) != 1) {
 489         if (i || m)
 490           fprintf(stderr,"error reading statusfile: x[%d][%d]\n",i,m);
 491         /* else eof! */
 492         return FALSE;
 493       }
 494       x[i][m]=x0;
 495     }
 496   }
 497   if (bReadBox) {
 498     for(m=0; (m<DIM); m++) {
 499       if (fscanf(status,"%lf",&x0) != 1)
 500         return FALSE;
 501       box[m][m]=x0;
 502     }
 503   }
 504   return TRUE;
 505 }
 506
 507 static bool xyz_next_x(FILE *status, const output_env_t oenv,
 508                        real *t, int natoms, rvec x[], matrix box)
 509      /* Reads until a new x can be found (return TRUE)
 510       * or eof (return FALSE)
 511       */
 512 {
 513   real pt;
 514
 515   pt=*t;
 516   while (!bTimeSet(TBEGIN) || (*t < rTimeValue(TBEGIN))) {
 517     if (!do_read_xyz(status,natoms,x,box))
 518       return FALSE;
 519     printcount(oenv,*t,FALSE);
 520     *t+=DT;
 521     pt=*t;
 522   }
 523   if (!bTimeSet(TEND) || (*t <= rTimeValue(TEND))) {
 524     if (!do_read_xyz(status,natoms,x,box)) {
 525       printlast(oenv,*t);
 526       return FALSE;
 527     }
 528     printcount(oenv,*t,FALSE);
 529     pt=*t;
 530     *t+=DT;
 531     return TRUE;
 532   }
 533   printlast(oenv,pt);
 534   return FALSE;
 535 }
 536
 537 static int xyz_first_x(FILE *status, const output_env_t oenv,
 538                        real *t, rvec **x, matrix box)
 539 /* Reads status, mallocs x, and returns x and box
 540  * Returns natoms when succesful, FALSE otherwise
 541  */
 542 {
 543   int    m;
 544
 545   INITCOUNT;
 546
 547   clear_mat(box);
 548   choose_ff(status);
 549
 550   for(m=0; (m<DIM); m++)
 551     box[m][m]=BOX[m];
 552
 553   snew(*x,NATOMS);
 554   *t=DT;
 555   if (!xyz_next_x(status,oenv,t,NATOMS,*x,box))
 556     return 0;
 557   *t=0.0;
 558
 559   return NATOMS;
 560 }
 561
 562 static bool pdb_next_x(FILE *status,t_trxframe *fr)
 563 {
 564   t_atoms   atoms;
 565   matrix    boxpdb;
 566   int       ePBC,model_nr,na;
 567   char      title[STRLEN],*time;
 568   double    dbl;
 569
 570   atoms.nr = fr->natoms;
 571   atoms.atom=NULL;
 572   atoms.pdbinfo=NULL;
 573   /* the other pointers in atoms should not be accessed if these are NULL */
 574   model_nr=NOTSET;
 575   na=read_pdbfile(status,title,&model_nr,&atoms,fr->x,&ePBC,boxpdb,TRUE,NULL);
 576   set_trxframe_ePBC(fr,ePBC);
 577   if (nframes_read()==0)
 578     fprintf(stderr," '%s', %d atoms\n",title, fr->natoms);
 579   fr->bPrec = TRUE;
 580   fr->prec = 10000;
 581   fr->bX = TRUE;
 582   fr->bBox = (boxpdb[XX][XX] != 0.0);
 583   if (fr->bBox) {
 584     copy_mat(boxpdb,fr->box);
 585   }
 586
 587   if (model_nr!=NOTSET) {
 588     fr->bStep = TRUE;
 589     fr->step = model_nr;
 590   }
 591   time=strstr(title," t= ");
 592   if (time) {
 593     fr->bTime = TRUE;
 594     sscanf(time+4,"%lf",&dbl);
 595     fr->time=(real)dbl;
 596   } else {
 597     fr->bTime = FALSE;
 598     /* this is a bit dirty, but it will work: if no time is read from
 599        comment line in pdb file, set time to current frame number */
 600     if (fr->bStep)
 601       fr->time=(real)fr->step;
 602     else
 603       fr->time=(real)nframes_read();
 604   }
 605   if (na == 0) {
 606     return FALSE;
 607   } else {
 608     if (na != fr->natoms)
 609       gmx_fatal(FARGS,"Number of atoms in pdb frame %d is %d instead of %d",
 610                   nframes_read(),na,fr->natoms);
 611     return TRUE;
 612   }
 613 }
 614
 615 static int pdb_first_x(FILE *status, t_trxframe *fr)
 616 {
 617   INITCOUNT;
 618
 619   fprintf(stderr,"Reading frames from pdb file");
 620   frewind(status);
 621   get_pdb_coordnum(status, &fr->natoms);
 622   if (fr->natoms==0)
 623     gmx_fatal(FARGS,"\nNo coordinates in pdb file\n");
 624   frewind(status);
 625   snew(fr->x,fr->natoms);
 626   pdb_next_x(status, fr);
 627
 628   return fr->natoms;
 629 }
 630
 631 bool read_next_frame(const output_env_t oenv,int status,t_trxframe *fr)
 632 {
 633   real pt;
 634   int  ct;
 635   bool bOK,bRet,bMissingData=FALSE,bSkip=FALSE;
 636
 637   bRet = FALSE;
 638   pt=fr->time;
 639
 640   do {
 641     clear_trxframe(fr,FALSE);
 642     fr->tppf = fr->tpf;
 643     fr->tpf  = fr->time;
 644
 645     switch (gmx_fio_getftp(status)) {
 646     case efTRJ:
 647     case efTRR:
 648         bRet = gmx_next_frame(status,fr);
 649         break;
 650     case efCPT:
 651       /* Checkpoint files can not contain mulitple frames */
 652       break;
 653     case efG96:
 654       read_g96_conf(gmx_fio_getfp(status),NULL,fr);
 655       bRet = (fr->natoms > 0);
 656       break;
 657     case efG87:
 658       bRet = xyz_next_x(gmx_fio_getfp(status),oenv,&fr->time,fr->natoms,
 659                         fr->x,fr->box);
 660       fr->bTime = bRet;
 661       fr->bX    = bRet;
 662       fr->bBox  = bRet;
 663       break;
 664     case efXTC:
 665       /* B. Hess 2005-4-20
 666        * Sometimes is off by one frame
 667        * and sometimes reports frame not present/file not seekable
 668        */
 669       /* DvdS 2005-05-31: this has been fixed along with the increased
 670        * accuracy of the control over -b and -e options.
 671        */
 672         if (bTimeSet(TBEGIN) && (fr->time < rTimeValue(TBEGIN))) {
 673             if (xtc_seek_time(rTimeValue(TBEGIN),status,fr->natoms)) {
 674                 gmx_fatal(FARGS,"Specified frame doesn't exist or file not seekable");
 675             }
 676             INITCOUNT;
 677         }
 678       bRet = read_next_xtc(status,fr->natoms,&fr->step,&fr->time,fr->box,
 679                            fr->x,&fr->prec,&bOK);
 680       fr->bPrec = (bRet && fr->prec > 0);
 681       fr->bStep = bRet;
 682       fr->bTime = bRet;
 683       fr->bX    = bRet;
 684       fr->bBox  = bRet;
 685       if (!bOK) {
 686         /* Actually the header could also be not ok,
 687            but from bOK from read_next_xtc this can't be distinguished */
 688         fr->not_ok = DATA_NOT_OK;
 689       }
 690       break;
 691     case efPDB:
 692       bRet = pdb_next_x(gmx_fio_getfp(status),fr);
 693       break;
 694     case efGRO:
 695       bRet = gro_next_x_or_v(gmx_fio_getfp(status),fr);
 696       break;
 697     default:
 698       gmx_fatal(FARGS,"DEATH HORROR in read_next_frame ftp=%s,status=%d",
 699                   ftp2ext(gmx_fio_getftp(status)),status);
 700     }
 701
 702     if (bRet) {
 703       bMissingData = ((fr->flags & TRX_NEED_X && !fr->bX) ||
 704                       (fr->flags & TRX_NEED_V && !fr->bV) ||
 705                       (fr->flags & TRX_NEED_F && !fr->bF));
 706       bSkip = FALSE;
 707       if (!bMissingData) {
 708         ct=check_times2(fr->time,fr->t0,fr->tpf,fr->tppf,fr->bDouble);
 709         if (ct == 0 || (fr->flags & TRX_DONT_SKIP && ct<0)) {
 710           printcount(oenv,fr->time,FALSE);
 711         } else if (ct > 0)
 712           bRet = FALSE;
 713         else {
 714           printcount(oenv,fr->time,TRUE);
 715           bSkip = TRUE;
 716         }
 717       }
 718     }
 719
 720   } while (bRet && (bMissingData || bSkip));
 721
 722   if (!bRet) {
 723     printlast(oenv,pt);
 724     if (fr->not_ok)
 725       printincomp(fr);
 726   }
 727
 728   return bRet;
 729 }
 730
 731 int read_first_frame(const output_env_t oenv,int *status,
 732                      const char *fn,t_trxframe *fr,int flags)
 733 {
 734   int  fp;
 735   bool bFirst,bOK;
 736
 737   clear_trxframe(fr,TRUE);
 738   fr->flags = flags;
 739
 740   bFirst = TRUE;
 741   INITCOUNT;
 742
 743   fp = *status =gmx_fio_open(fn,"r");
 744   switch (gmx_fio_getftp(fp))
 745   {
 746   case efTRJ:
 747   case efTRR:
 748     break;
 749   case efCPT:
 750     read_checkpoint_trxframe(fp,fr);
 751     bFirst = FALSE;
 752     break;
 753   case efG96:
 754     /* Can not rewind a compressed file, so open it twice */
 755     read_g96_conf(gmx_fio_getfp(fp),fn,fr);
 756     gmx_fio_close(fp);
 757     clear_trxframe(fr,FALSE);
 758     if (flags & (TRX_READ_X | TRX_NEED_X))
 759       snew(fr->x,fr->natoms);
 760     if (flags & (TRX_READ_V | TRX_NEED_V))
 761       snew(fr->v,fr->natoms);
 762     fp = *status =gmx_fio_open(fn,"r");
 763     break;
 764   case efG87:
 765     fr->natoms=xyz_first_x(gmx_fio_getfp(fp),oenv,&fr->time,&fr->x,fr->box);
 766     if (fr->natoms) {
 767       fr->bTime = TRUE;
 768       fr->bX    = TRUE;
 769       fr->bBox  = TRUE;
 770       printcount(oenv,fr->time,FALSE);
 771     }
 772     bFirst = FALSE;
 773     break;
 774   case efXTC:
 775     if (read_first_xtc(fp,&fr->natoms,&fr->step,&fr->time,fr->box,&fr->x,
 776                        &fr->prec,&bOK) == 0) {
 777       if (bOK) {
 778         gmx_fatal(FARGS,"No XTC!\n");
 779       } else {
 780         fr->not_ok = DATA_NOT_OK;
 781       }
 782     }
 783     if (fr->not_ok) {
 784       fr->natoms = 0;
 785       printincomp(fr);
 786     } else {
 787       fr->bPrec = (fr->prec > 0);
 788       fr->bStep = TRUE;
 789       fr->bTime = TRUE;
 790       fr->bX    = TRUE;
 791       fr->bBox  = TRUE;
 792       printcount(oenv,fr->time,FALSE);
 793     }
 794     bFirst = FALSE;
 795     break;
 796   case efPDB:
 797     pdb_first_x(gmx_fio_getfp(fp),fr);
 798     if (fr->natoms)
 799       printcount(oenv,fr->time,FALSE);
 800     bFirst = FALSE;
 801     break;
 802   case efGRO:
 803     if (gro_first_x_or_v(gmx_fio_getfp(fp),fr))
 804       printcount(oenv,fr->time,FALSE);
 805     bFirst = FALSE;
 806     break;
 807   default:
 808     gmx_fatal(FARGS,"Not supported in read_first_frame: %s",fn);
 809     break;
 810   }
 811
 812   if (bFirst ||
 813       (!(fr->flags & TRX_DONT_SKIP) && check_times(fr->time) < 0))
 814     /* Read a frame when no frame was read or the first was skipped */
 815     if (!read_next_frame(oenv,*status,fr))
 816       return FALSE;
 817   fr->t0 = fr->time;
 818
 819   return (fr->natoms > 0);
 820 }
 821
 822 /***** C O O R D I N A T E   S T U F F *****/
 823
 824 int read_first_x(const output_env_t oenv,int *status,const char *fn,
 825                  real *t,rvec **x,matrix box)
 826 {
 827   t_trxframe fr;
 828
 829   read_first_frame(oenv,status,fn,&fr,TRX_NEED_X);
 830   if (*status >= nxframe) {
 831     nxframe = *status+1;
 832     srenew(xframe,nxframe);
 833   }
 834   xframe[*status] = fr;
 835   *t = xframe[*status].time;
 836   *x = xframe[*status].x;
 837   copy_mat(xframe[*status].box,box);
 838
 839   return xframe[*status].natoms;
 840 }
 841
 842 bool read_next_x(const output_env_t oenv, int status,real *t, int natoms,
 843                  rvec x[], matrix box)
 844 {
 845   bool bRet;
 846
 847   xframe[status].x = x;
 848   bRet = read_next_frame(oenv,status,&xframe[status]);
 849   *t = xframe[status].time;
 850   copy_mat(xframe[status].box,box);
 851
 852   return bRet;
 853 }
 854
 855 void close_trj(int status)
 856 {
 857   gmx_fio_close(status);
 858 }
 859
 860 void rewind_trj(int status)
 861 {
 862   INITCOUNT;
 863
 864   gmx_fio_rewind(status);
 865 }
 866
 867 /***** V E L O C I T Y   S T U F F *****/
 868
 869 static void clear_v(t_trxframe *fr)
 870 {
 871   int i;
 872
 873   if (!fr->bV)
 874     for(i=0; i<fr->natoms; i++)
 875       clear_rvec(fr->v[i]);
 876 }
 877
 878 int read_first_v(const output_env_t oenv, int *status,const char *fn,real *t,
 879                  rvec **v,matrix box)
 880 {
 881   t_trxframe fr;
 882
 883   read_first_frame(oenv,status,fn,&fr,TRX_NEED_V);
 884   *t = fr.time;
 885   clear_v(&fr);
 886   *v = fr.v;
 887   copy_mat(fr.box,box);
 888
 889   return fr.natoms;
 890 }
 891
 892 bool read_next_v(const output_env_t oenv,int status,real *t,int natoms,rvec v[],
 893                  matrix box)
 894 {
 895   t_trxframe fr;
 896   bool bRet;
 897
 898   clear_trxframe(&fr,TRUE);
 899   fr.flags = TRX_NEED_V;
 900   fr.natoms = natoms;
 901   fr.time = *t;
 902   fr.v = v;
 903   bRet = read_next_frame(oenv,status,&fr);
 904   *t = fr.time;
 905   clear_v(&fr);
 906   copy_mat(fr.box,box);
 907
 908   return bRet;
 909 }