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