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