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