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