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