src/gmxlib/trxio.c

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