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