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