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