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