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     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         output_env_get_trajectory_io_verbosity(oenv) != 0)
 193     {
 194         fprintf(stderr, "\r%-14s %6d time %8.3f   ", l, status->__frame,
 195                 output_env_conv_time(oenv, t));
 196         fflush(stderr);
 197     }
 198 }
 199
 200 static void printcount(t_trxstatus *status, const gmx_output_env_t *oenv, real t,
 201                        gmx_bool bSkip)
 202 {
 203     status->__frame++;
 204     printcount_(status, oenv, bSkip ? "Skipping frame" : "Reading frame", t);
 205 }
 206
 207 static void printlast(t_trxstatus *status, const gmx_output_env_t *oenv, real t)
 208 {
 209     printcount_(status, oenv, "Last frame", t);
 210     fprintf(stderr, "\n");
 211     fflush(stderr);
 212 }
 213
 214 static void printincomp(t_trxstatus *status, t_trxframe *fr)
 215 {
 216     if (fr->not_ok & HEADER_NOT_OK)
 217     {
 218         fprintf(stderr, "WARNING: Incomplete header: nr %d time %g\n",
 219                 status->__frame+1, fr->time);
 220     }
 221     else if (fr->not_ok)
 222     {
 223         fprintf(stderr, "WARNING: Incomplete frame: nr %d time %g\n",
 224                 status->__frame+1, fr->time);
 225     }
 226     fflush(stderr);
 227 }
 228
 229 int prec2ndec(real prec)
 230 {
 231     if (prec <= 0)
 232     {
 233         gmx_fatal(FARGS, "DEATH HORROR prec (%g) <= 0 in prec2ndec", prec);
 234     }
 235
 236     return (int)(log(prec)/log(10.0)+0.5);
 237 }
 238
 239 real ndec2prec(int ndec)
 240 {
 241     return pow(10.0, ndec);
 242 }
 243
 244 t_fileio *trx_get_fileio(t_trxstatus *status)
 245 {
 246     return status->fio;
 247 }
 248
 249 float trx_get_time_of_final_frame(t_trxstatus *status)
 250 {
 251     t_fileio *stfio    = trx_get_fileio(status);
 252     int       filetype = gmx_fio_getftp(stfio);
 253     int       bOK;
 254     float     lasttime = -1;
 255
 256     if (filetype == efXTC)
 257     {
 258         lasttime =
 259             xdr_xtc_get_last_frame_time(gmx_fio_getfp(stfio),
 260                                         gmx_fio_getxdr(stfio),
 261                                         status->natoms, &bOK);
 262         if (!bOK)
 263         {
 264             gmx_fatal(FARGS, "Error reading last frame. Maybe seek not supported." );
 265         }
 266     }
 267     else if (filetype == efTNG)
 268     {
 269         gmx_tng_trajectory_t tng = status->tng;
 270         if (!tng)
 271         {
 272             gmx_fatal(FARGS, "Error opening TNG file.");
 273         }
 274         lasttime = gmx_tng_get_time_of_final_frame(tng);
 275     }
 276     else
 277     {
 278         gmx_incons("Only supported for TNG and XTC");
 279     }
 280     return lasttime;
 281 }
 282
 283 void clear_trxframe(t_trxframe *fr, gmx_bool bFirst)
 284 {
 285     fr->not_ok    = 0;
 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->step      = 0;
 301         fr->time      = 0;
 302         fr->lambda    = 0;
 303         fr->fep_state = 0;
 304         fr->atoms     = nullptr;
 305         fr->prec      = 0;
 306         fr->x         = nullptr;
 307         fr->v         = nullptr;
 308         fr->f         = nullptr;
 309         clear_mat(fr->box);
 310         fr->bPBC   = FALSE;
 311         fr->ePBC   = -1;
 312     }
 313 }
 314
 315 void set_trxframe_ePBC(t_trxframe *fr, int ePBC)
 316 {
 317     fr->bPBC = (ePBC == -1);
 318     fr->ePBC = ePBC;
 319 }
 320
 321 int write_trxframe_indexed(t_trxstatus *status, const t_trxframe *fr, int nind,
 322                            const int *ind, gmx_conect gc)
 323 {
 324     char  title[STRLEN];
 325     rvec *xout = nullptr, *vout = nullptr, *fout = nullptr;
 326     int   i, ftp = -1;
 327     real  prec;
 328
 329     if (fr->bPrec)
 330     {
 331         prec = fr->prec;
 332     }
 333     else
 334     {
 335         prec = 1000.0;
 336     }
 337
 338     if (status->tng)
 339     {
 340         ftp = efTNG;
 341     }
 342     else if (status->fio)
 343     {
 344         ftp = gmx_fio_getftp(status->fio);
 345     }
 346     else
 347     {
 348         gmx_incons("No input file available");
 349     }
 350
 351     switch (ftp)
 352     {
 353         case efTRR:
 354         case efTNG:
 355             break;
 356         default:
 357             if (!fr->bX)
 358             {
 359                 gmx_fatal(FARGS, "Need coordinates to write a %s trajectory",
 360                           ftp2ext(ftp));
 361             }
 362             break;
 363     }
 364
 365     switch (ftp)
 366     {
 367         case efTRR:
 368         case efTNG:
 369             if (fr->bV)
 370             {
 371                 snew(vout, nind);
 372                 for (i = 0; i < nind; i++)
 373                 {
 374                     copy_rvec(fr->v[ind[i]], vout[i]);
 375                 }
 376             }
 377             if (fr->bF)
 378             {
 379                 snew(fout, nind);
 380                 for (i = 0; i < nind; i++)
 381                 {
 382                     copy_rvec(fr->f[ind[i]], fout[i]);
 383                 }
 384             }
 385         // fallthrough
 386         case efXTC:
 387             if (fr->bX)
 388             {
 389                 snew(xout, nind);
 390                 for (i = 0; i < nind; i++)
 391                 {
 392                     copy_rvec(fr->x[ind[i]], xout[i]);
 393                 }
 394             }
 395             break;
 396         default:
 397             break;
 398     }
 399
 400     switch (ftp)
 401     {
 402         case efTNG:
 403             gmx_write_tng_from_trxframe(status->tng, fr, nind);
 404             break;
 405         case efXTC:
 406             write_xtc(status->fio, nind, fr->step, fr->time, fr->box, xout, prec);
 407             break;
 408         case efTRR:
 409             gmx_trr_write_frame(status->fio, nframes_read(status),
 410                                 fr->time, fr->step, fr->box, nind, xout, vout, fout);
 411             break;
 412         case efGRO:
 413         case efPDB:
 414         case efBRK:
 415         case efENT:
 416             if (!fr->bAtoms)
 417             {
 418                 gmx_fatal(FARGS, "Can not write a %s file without atom names",
 419                           ftp2ext(ftp));
 420             }
 421             sprintf(title, "frame t= %.3f", fr->time);
 422             if (ftp == efGRO)
 423             {
 424                 write_hconf_indexed_p(gmx_fio_getfp(status->fio), title, fr->atoms, nind, ind,
 425                                       fr->x, fr->bV ? fr->v : nullptr, fr->box);
 426             }
 427             else
 428             {
 429                 write_pdbfile_indexed(gmx_fio_getfp(status->fio), title, fr->atoms,
 430                                       fr->x, -1, fr->box, ' ', fr->step, nind, ind, gc, TRUE);
 431             }
 432             break;
 433         case efG96:
 434             sprintf(title, "frame t= %.3f", fr->time);
 435             write_g96_conf(gmx_fio_getfp(status->fio), title, 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         // fallthrough
 456         case efXTC:
 457             sfree(xout);
 458             break;
 459         default:
 460             break;
 461     }
 462
 463     return 0;
 464 }
 465
 466 t_trxstatus *
 467 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                                  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(filename,
 487                                 filemode,
 488                                 &in->tng,
 489                                 &out->tng,
 490                                 natoms,
 491                                 mtop,
 492                                 index,
 493                                 index_group_name);
 494     }
 495     else if ((infile) && (efTNG == fn2ftp(infile)))
 496     {
 497         gmx_tng_trajectory_t tng_in;
 498         gmx_tng_open(infile, 'r', &tng_in);
 499
 500         gmx_prepare_tng_writing(filename,
 501                                 filemode,
 502                                 &tng_in,
 503                                 &out->tng,
 504                                 natoms,
 505                                 mtop,
 506                                 index,
 507                                 index_group_name);
 508     }
 509     else
 510     {
 511         // we start from a file that is not a tng file or have been unable to load the
 512         // input file, so we need to populate the fields independently of it
 513         gmx_prepare_tng_writing(filename,
 514                                 filemode,
 515                                 nullptr,
 516                                 &out->tng,
 517                                 natoms,
 518                                 mtop,
 519                                 index,
 520                                 index_group_name);
 521     }
 522     return out;
 523 }
 524
 525 void write_tng_frame(t_trxstatus *status,
 526                      t_trxframe  *frame)
 527 {
 528     gmx_write_tng_from_trxframe(status->tng, frame, -1);
 529 }
 530
 531 int write_trxframe(t_trxstatus *status, t_trxframe *fr, gmx_conect gc)
 532 {
 533     char title[STRLEN];
 534     real prec;
 535
 536     if (fr->bPrec)
 537     {
 538         prec = fr->prec;
 539     }
 540     else
 541     {
 542         prec = 1000.0;
 543     }
 544
 545     if (status->tng)
 546     {
 547         gmx_tng_set_compression_precision(status->tng, prec);
 548         write_tng_frame(status, fr);
 549
 550         return 0;
 551     }
 552
 553     switch (gmx_fio_getftp(status->fio))
 554     {
 555         case efTRR:
 556             break;
 557         default:
 558             if (!fr->bX)
 559             {
 560                 gmx_fatal(FARGS, "Need coordinates to write a %s trajectory",
 561                           ftp2ext(gmx_fio_getftp(status->fio)));
 562             }
 563             break;
 564     }
 565
 566     switch (gmx_fio_getftp(status->fio))
 567     {
 568         case efXTC:
 569             write_xtc(status->fio, fr->natoms, fr->step, fr->time, fr->box, fr->x, prec);
 570             break;
 571         case efTRR:
 572             gmx_trr_write_frame(status->fio, fr->step, fr->time, fr->lambda, fr->box, fr->natoms,
 573                                 fr->bX ? fr->x : nullptr, fr->bV ? fr->v : nullptr, fr->bF ? fr->f : nullptr);
 574             break;
 575         case efGRO:
 576         case efPDB:
 577         case efBRK:
 578         case efENT:
 579             if (!fr->bAtoms)
 580             {
 581                 gmx_fatal(FARGS, "Can not write a %s file without atom names",
 582                           ftp2ext(gmx_fio_getftp(status->fio)));
 583             }
 584             sprintf(title, "frame t= %.3f", fr->time);
 585             if (gmx_fio_getftp(status->fio) == efGRO)
 586             {
 587                 write_hconf_p(gmx_fio_getfp(status->fio), title, fr->atoms,
 588                               fr->x, fr->bV ? fr->v : nullptr, fr->box);
 589             }
 590             else
 591             {
 592                 write_pdbfile(gmx_fio_getfp(status->fio), title,
 593                               fr->atoms, fr->x, fr->bPBC ? fr->ePBC : -1, fr->box,
 594                               ' ', fr->step, gc, TRUE);
 595             }
 596             break;
 597         case efG96:
 598             write_g96_conf(gmx_fio_getfp(status->fio), title, fr, -1, nullptr);
 599             break;
 600         default:
 601             gmx_fatal(FARGS, "Sorry, write_trxframe can not write %s",
 602                       ftp2ext(gmx_fio_getftp(status->fio)));
 603             break;
 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, "%" GMX_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 }