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