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