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