src/gromacs/fileio/trxio.cpp

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