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