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