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