src/gromacs/fileio/trxio.c

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