src/gromacs/fileio/trxio.cpp

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