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