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 "trxio.h"
  38
  39 #ifdef HAVE_CONFIG_H
  40 #include <config.h>
  41 #endif
  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 "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 "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->tpf       = 0;
 285         fr->tppf      = 0;
 286         fr->title     = NULL;
 287         fr->step      = 0;
 288         fr->time      = 0;
 289         fr->lambda    = 0;
 290         fr->fep_state = 0;
 291         fr->atoms     = NULL;
 292         fr->prec      = 0;
 293         fr->x         = NULL;
 294         fr->v         = NULL;
 295         fr->f         = NULL;
 296         clear_mat(fr->box);
 297         fr->bPBC   = FALSE;
 298         fr->ePBC   = -1;
 299     }
 300 }
 301
 302 void set_trxframe_ePBC(t_trxframe *fr, int ePBC)
 303 {
 304     fr->bPBC = (ePBC == -1);
 305     fr->ePBC = ePBC;
 306 }
 307
 308 int write_trxframe_indexed(t_trxstatus *status, t_trxframe *fr, int nind,
 309                            const atom_id *ind, gmx_conect gc)
 310 {
 311     char  title[STRLEN];
 312     rvec *xout = NULL, *vout = NULL, *fout = NULL;
 313     int   i, ftp = -1;
 314     real  prec;
 315
 316     if (fr->bPrec)
 317     {
 318         prec = fr->prec;
 319     }
 320     else
 321     {
 322         prec = 1000.0;
 323     }
 324
 325     if (status->tng)
 326     {
 327         ftp = efTNG;
 328     }
 329     else if (status->fio)
 330     {
 331         ftp = gmx_fio_getftp(status->fio);
 332     }
 333     else
 334     {
 335         gmx_incons("No input file available");
 336     }
 337
 338     switch (ftp)
 339     {
 340         case efTRJ:
 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 efTRJ:
 356         case efTRR:
 357         case efTNG:
 358             if (fr->bV)
 359             {
 360                 snew(vout, nind);
 361                 for (i = 0; i < nind; i++)
 362                 {
 363                     copy_rvec(fr->v[ind[i]], vout[i]);
 364                 }
 365             }
 366             if (fr->bF)
 367             {
 368                 snew(fout, nind);
 369                 for (i = 0; i < nind; i++)
 370                 {
 371                     copy_rvec(fr->f[ind[i]], fout[i]);
 372                 }
 373             }
 374         /* no break */
 375         case efXTC:
 376             if (fr->bX)
 377             {
 378                 snew(xout, nind);
 379                 for (i = 0; i < nind; i++)
 380                 {
 381                     copy_rvec(fr->x[ind[i]], xout[i]);
 382                 }
 383             }
 384             break;
 385         default:
 386             break;
 387     }
 388
 389     switch (ftp)
 390     {
 391         case efTNG:
 392             gmx_write_tng_from_trxframe(status->tng, fr, nind);
 393             break;
 394         case efXTC:
 395             write_xtc(status->fio, nind, fr->step, fr->time, fr->box, xout, prec);
 396             break;
 397         case efTRJ:
 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 efTRN:
 436         case efTRJ:
 437         case efTRR:
 438         case efTNG:
 439             if (vout)
 440             {
 441                 sfree(vout);
 442             }
 443             if (fout)
 444             {
 445                 sfree(fout);
 446             }
 447         /* no break */
 448         case efXTC:
 449             sfree(xout);
 450             break;
 451         default:
 452             break;
 453     }
 454
 455     return 0;
 456 }
 457
 458 void trjtools_gmx_prepare_tng_writing(const char       *filename,
 459                                       char              filemode,
 460                                       t_trxstatus      *in,
 461                                       t_trxstatus     **out,
 462                                       const char       *infile,
 463                                       const int         natoms,
 464                                       const gmx_mtop_t *mtop,
 465                                       const atom_id    *index,
 466                                       const char       *index_group_name)
 467 {
 468     if (filemode != 'w' && filemode != 'a')
 469     {
 470         gmx_incons("Sorry, can only prepare for TNG output.");
 471     }
 472
 473     if (*out == NULL)
 474     {
 475         snew((*out), 1);
 476     }
 477     status_init(*out);
 478
 479     if (in != NULL)
 480     {
 481         gmx_prepare_tng_writing(filename,
 482                                 filemode,
 483                                 &in->tng,
 484                                 &(*out)->tng,
 485                                 natoms,
 486                                 mtop,
 487                                 index,
 488                                 index_group_name);
 489     }
 490     else if (efTNG == fn2ftp(infile))
 491     {
 492         tng_trajectory_t tng_in;
 493         gmx_tng_open(infile, 'r', &tng_in);
 494
 495         gmx_prepare_tng_writing(filename,
 496                                 filemode,
 497                                 &tng_in,
 498                                 &(*out)->tng,
 499                                 natoms,
 500                                 mtop,
 501                                 index,
 502                                 index_group_name);
 503     }
 504 }
 505
 506 void write_tng_frame(t_trxstatus *status,
 507                      t_trxframe  *frame)
 508 {
 509     gmx_write_tng_from_trxframe(status->tng, frame, -1);
 510 }
 511
 512 int write_trxframe(t_trxstatus *status, t_trxframe *fr, gmx_conect gc)
 513 {
 514     char title[STRLEN];
 515     real prec;
 516
 517     if (fr->bPrec)
 518     {
 519         prec = fr->prec;
 520     }
 521     else
 522     {
 523         prec = 1000.0;
 524     }
 525
 526     if (status->tng)
 527     {
 528         gmx_tng_set_compression_precision(status->tng, prec);
 529         write_tng_frame(status, fr);
 530
 531         return 0;
 532     }
 533
 534     switch (gmx_fio_getftp(status->fio))
 535     {
 536         case efTRJ:
 537         case efTRR:
 538             break;
 539         default:
 540             if (!fr->bX)
 541             {
 542                 gmx_fatal(FARGS, "Need coordinates to write a %s trajectory",
 543                           ftp2ext(gmx_fio_getftp(status->fio)));
 544             }
 545             break;
 546     }
 547
 548     switch (gmx_fio_getftp(status->fio))
 549     {
 550         case efXTC:
 551             write_xtc(status->fio, fr->natoms, fr->step, fr->time, fr->box, fr->x, prec);
 552             break;
 553         case efTRJ:
 554         case efTRR:
 555             fwrite_trn(status->fio, fr->step, fr->time, fr->lambda, fr->box, fr->natoms,
 556                        fr->bX ? fr->x : NULL, fr->bV ? fr->v : NULL, fr->bF ? fr->f : NULL);
 557             break;
 558         case efGRO:
 559         case efPDB:
 560         case efBRK:
 561         case efENT:
 562             if (!fr->bAtoms)
 563             {
 564                 gmx_fatal(FARGS, "Can not write a %s file without atom names",
 565                           ftp2ext(gmx_fio_getftp(status->fio)));
 566             }
 567             sprintf(title, "frame t= %.3f", fr->time);
 568             if (gmx_fio_getftp(status->fio) == efGRO)
 569             {
 570                 write_hconf_p(gmx_fio_getfp(status->fio), title, fr->atoms,
 571                               prec2ndec(prec), fr->x, fr->bV ? fr->v : NULL, fr->box);
 572             }
 573             else
 574             {
 575                 write_pdbfile(gmx_fio_getfp(status->fio), title,
 576                               fr->atoms, fr->x, fr->bPBC ? fr->ePBC : -1, fr->box,
 577                               ' ', fr->step, gc, TRUE);
 578             }
 579             break;
 580         case efG96:
 581             write_g96_conf(gmx_fio_getfp(status->fio), fr, -1, NULL);
 582             break;
 583         default:
 584             gmx_fatal(FARGS, "Sorry, write_trxframe can not write %s",
 585                       ftp2ext(gmx_fio_getftp(status->fio)));
 586             break;
 587     }
 588
 589     return 0;
 590 }
 591
 592 int write_trx(t_trxstatus *status, int nind, const atom_id *ind, t_atoms *atoms,
 593               int step, real time, matrix box, rvec x[], rvec *v,
 594               gmx_conect gc)
 595 {
 596     t_trxframe fr;
 597
 598     clear_trxframe(&fr, TRUE);
 599     fr.bStep  = TRUE;
 600     fr.step   = step;
 601     fr.bTime  = TRUE;
 602     fr.time   = time;
 603     fr.bAtoms = atoms != NULL;
 604     fr.atoms  = atoms;
 605     fr.bX     = TRUE;
 606     fr.x      = x;
 607     fr.bV     = v != NULL;
 608     fr.v      = v;
 609     fr.bBox   = TRUE;
 610     copy_mat(box, fr.box);
 611
 612     return write_trxframe_indexed(status, &fr, nind, ind, gc);
 613 }
 614
 615 void close_trx(t_trxstatus *status)
 616 {
 617     gmx_tng_close(&status->tng);
 618     if (status->fio)
 619     {
 620         gmx_fio_close(status->fio);
 621     }
 622     sfree(status);
 623 }
 624
 625 t_trxstatus *open_trx(const char *outfile, const char *filemode)
 626 {
 627     t_trxstatus *stat;
 628     if (filemode[0] != 'w' && filemode[0] != 'a' && filemode[1] != '+')
 629     {
 630         gmx_fatal(FARGS, "Sorry, write_trx can only write");
 631     }
 632
 633     snew(stat, 1);
 634     status_init(stat);
 635
 636     stat->fio = gmx_fio_open(outfile, filemode);
 637     return stat;
 638 }
 639
 640 static gmx_bool gmx_next_frame(t_trxstatus *status, t_trxframe *fr)
 641 {
 642     t_trnheader sh;
 643     gmx_bool    bOK, bRet;
 644
 645     bRet = FALSE;
 646
 647     if (fread_trnheader(status->fio, &sh, &bOK))
 648     {
 649         fr->bDouble   = sh.bDouble;
 650         fr->natoms    = sh.natoms;
 651         fr->bStep     = TRUE;
 652         fr->step      = sh.step;
 653         fr->bTime     = TRUE;
 654         fr->time      = sh.t;
 655         fr->bLambda   = TRUE;
 656         fr->bFepState = TRUE;
 657         fr->lambda    = sh.lambda;
 658         fr->bBox      = sh.box_size > 0;
 659         if (fr->flags & (TRX_READ_X | TRX_NEED_X))
 660         {
 661             if (fr->x == NULL)
 662             {
 663                 snew(fr->x, sh.natoms);
 664             }
 665             fr->bX = sh.x_size > 0;
 666         }
 667         if (fr->flags & (TRX_READ_V | TRX_NEED_V))
 668         {
 669             if (fr->v == NULL)
 670             {
 671                 snew(fr->v, sh.natoms);
 672             }
 673             fr->bV = sh.v_size > 0;
 674         }
 675         if (fr->flags & (TRX_READ_F | TRX_NEED_F))
 676         {
 677             if (fr->f == NULL)
 678             {
 679                 snew(fr->f, sh.natoms);
 680             }
 681             fr->bF = sh.f_size > 0;
 682         }
 683         if (fread_htrn(status->fio, &sh, fr->box, fr->x, fr->v, fr->f))
 684         {
 685             bRet = TRUE;
 686         }
 687         else
 688         {
 689             fr->not_ok = DATA_NOT_OK;
 690         }
 691     }
 692     else
 693     if (!bOK)
 694     {
 695         fr->not_ok = HEADER_NOT_OK;
 696     }
 697
 698     return bRet;
 699 }
 700
 701 static gmx_bool pdb_next_x(t_trxstatus *status, FILE *fp, t_trxframe *fr)
 702 {
 703     t_atoms   atoms;
 704     matrix    boxpdb;
 705     int       ePBC, model_nr, na;
 706     char      title[STRLEN], *time;
 707     double    dbl;
 708
 709     atoms.nr      = fr->natoms;
 710     atoms.atom    = NULL;
 711     atoms.pdbinfo = NULL;
 712     /* the other pointers in atoms should not be accessed if these are NULL */
 713     model_nr = NOTSET;
 714     na       = read_pdbfile(fp, title, &model_nr, &atoms, fr->x, &ePBC, boxpdb, TRUE, NULL);
 715     set_trxframe_ePBC(fr, ePBC);
 716     if (nframes_read(status) == 0)
 717     {
 718         fprintf(stderr, " '%s', %d atoms\n", title, fr->natoms);
 719     }
 720     fr->bPrec = TRUE;
 721     fr->prec  = 10000;
 722     fr->bX    = TRUE;
 723     fr->bBox  = (boxpdb[XX][XX] != 0.0);
 724     if (fr->bBox)
 725     {
 726         copy_mat(boxpdb, fr->box);
 727     }
 728
 729     if (model_nr != NOTSET)
 730     {
 731         fr->bStep = TRUE;
 732         fr->step  = model_nr;
 733     }
 734     time = strstr(title, " t= ");
 735     if (time)
 736     {
 737         fr->bTime = TRUE;
 738         sscanf(time+4, "%lf", &dbl);
 739         fr->time = (real)dbl;
 740     }
 741     else
 742     {
 743         fr->bTime = FALSE;
 744         /* this is a bit dirty, but it will work: if no time is read from
 745            comment line in pdb file, set time to current frame number */
 746         if (fr->bStep)
 747         {
 748             fr->time = (real)fr->step;
 749         }
 750         else
 751         {
 752             fr->time = (real)nframes_read(status);
 753         }
 754     }
 755     if (na == 0)
 756     {
 757         return FALSE;
 758     }
 759     else
 760     {
 761         if (na != fr->natoms)
 762         {
 763             gmx_fatal(FARGS, "Number of atoms in pdb frame %d is %d instead of %d",
 764                       nframes_read(status), na, fr->natoms);
 765         }
 766         return TRUE;
 767     }
 768 }
 769
 770 static int pdb_first_x(t_trxstatus *status, FILE *fp, t_trxframe *fr)
 771 {
 772     initcount(status);
 773
 774     fprintf(stderr, "Reading frames from pdb file");
 775     frewind(fp);
 776     get_pdb_coordnum(fp, &fr->natoms);
 777     if (fr->natoms == 0)
 778     {
 779         gmx_fatal(FARGS, "\nNo coordinates in pdb file\n");
 780     }
 781     frewind(fp);
 782     snew(fr->x, fr->natoms);
 783     pdb_next_x(status, fp, fr);
 784
 785     return fr->natoms;
 786 }
 787
 788 gmx_bool read_next_frame(const output_env_t oenv, t_trxstatus *status, t_trxframe *fr)
 789 {
 790     real     pt;
 791     int      ct;
 792     gmx_bool bOK, bRet, bMissingData = FALSE, bSkip = FALSE;
 793     int      dummy = 0;
 794     int      ftp;
 795
 796     bRet = FALSE;
 797     pt   = fr->time;
 798
 799     do
 800     {
 801         clear_trxframe(fr, FALSE);
 802         fr->tppf = fr->tpf;
 803         fr->tpf  = fr->time;
 804
 805         if (status->tng)
 806         {
 807             /* Special treatment for TNG files */
 808             ftp = efTNG;
 809         }
 810         else
 811         {
 812             ftp = gmx_fio_getftp(status->fio);
 813         }
 814         switch (ftp)
 815         {
 816             case efTRJ:
 817             case efTRR:
 818                 bRet = gmx_next_frame(status, fr);
 819                 break;
 820             case efCPT:
 821                 /* Checkpoint files can not contain mulitple frames */
 822                 break;
 823             case efG96:
 824                 read_g96_conf(gmx_fio_getfp(status->fio), NULL, fr,
 825                               status->persistent_line);
 826                 bRet = (fr->natoms > 0);
 827                 break;
 828             case efXTC:
 829                 /* B. Hess 2005-4-20
 830                  * Sometimes is off by one frame
 831                  * and sometimes reports frame not present/file not seekable
 832                  */
 833                 /* DvdS 2005-05-31: this has been fixed along with the increased
 834                  * accuracy of the control over -b and -e options.
 835                  */
 836                 if (bTimeSet(TBEGIN) && (fr->time < rTimeValue(TBEGIN)))
 837                 {
 838                     if (xtc_seek_time(status->fio, rTimeValue(TBEGIN), fr->natoms, TRUE))
 839                     {
 840                         gmx_fatal(FARGS, "Specified frame (time %f) doesn't exist or file corrupt/inconsistent.",
 841                                   rTimeValue(TBEGIN));
 842                     }
 843                     initcount(status);
 844                 }
 845                 bRet = read_next_xtc(status->fio, fr->natoms, &fr->step, &fr->time, fr->box,
 846                                      fr->x, &fr->prec, &bOK);
 847                 fr->bPrec = (bRet && fr->prec > 0);
 848                 fr->bStep = bRet;
 849                 fr->bTime = bRet;
 850                 fr->bX    = bRet;
 851                 fr->bBox  = bRet;
 852                 if (!bOK)
 853                 {
 854                     /* Actually the header could also be not ok,
 855                        but from bOK from read_next_xtc this can't be distinguished */
 856                     fr->not_ok = DATA_NOT_OK;
 857                 }
 858                 break;
 859             case efTNG:
 860                 bRet = gmx_read_next_tng_frame(status->tng, fr, NULL, 0);
 861                 break;
 862             case efPDB:
 863                 bRet = pdb_next_x(status, gmx_fio_getfp(status->fio), fr);
 864                 break;
 865             case efGRO:
 866                 bRet = gro_next_x_or_v(gmx_fio_getfp(status->fio), fr);
 867                 break;
 868             default:
 869 #ifdef GMX_USE_PLUGINS
 870                 bRet = read_next_vmd_frame(fr);
 871 #else
 872                 gmx_fatal(FARGS, "DEATH HORROR in read_next_frame ftp=%s,status=%s",
 873                           ftp2ext(gmx_fio_getftp(status->fio)),
 874                           gmx_fio_getname(status->fio));
 875 #endif
 876         }
 877
 878         if (bRet)
 879         {
 880             bMissingData = (((fr->flags & TRX_NEED_X) && !fr->bX) ||
 881                             ((fr->flags & TRX_NEED_V) && !fr->bV) ||
 882                             ((fr->flags & TRX_NEED_F) && !fr->bF));
 883             bSkip = FALSE;
 884             if (!bMissingData)
 885             {
 886                 ct = check_times2(fr->time, fr->t0, fr->bDouble);
 887                 if (ct == 0 || ((fr->flags & TRX_DONT_SKIP) && ct < 0))
 888                 {
 889                     printcount(status, oenv, fr->time, FALSE);
 890                 }
 891                 else if (ct > 0)
 892                 {
 893                     bRet = FALSE;
 894                 }
 895                 else
 896                 {
 897                     printcount(status, oenv, fr->time, TRUE);
 898                     bSkip = TRUE;
 899                 }
 900             }
 901         }
 902
 903     }
 904     while (bRet && (bMissingData || bSkip));
 905
 906     if (!bRet)
 907     {
 908         printlast(status, oenv, pt);
 909         if (fr->not_ok)
 910         {
 911             printincomp(status, fr);
 912         }
 913     }
 914
 915     return bRet;
 916 }
 917
 918 int read_first_frame(const output_env_t oenv, t_trxstatus **status,
 919                      const char *fn, t_trxframe *fr, int flags)
 920 {
 921     t_fileio      *fio;
 922     gmx_bool       bFirst, bOK;
 923     int            dummy = 0;
 924     int            ftp   = fn2ftp(fn);
 925     gmx_int64_t   *tng_ids;
 926
 927     clear_trxframe(fr, TRUE);
 928     fr->flags = flags;
 929
 930     bFirst = TRUE;
 931
 932     snew((*status), 1);
 933
 934     status_init( *status );
 935     (*status)->nxframe = 1;
 936     initcount(*status);
 937
 938     if (efTNG == ftp)
 939     {
 940         /* Special treatment for TNG files */
 941         gmx_tng_open(fn, 'r', &(*status)->tng);
 942     }
 943     else
 944     {
 945         fio = (*status)->fio = gmx_fio_open(fn, "r");
 946     }
 947     switch (ftp)
 948     {
 949         case efTRJ:
 950         case efTRR:
 951             break;
 952         case efCPT:
 953             read_checkpoint_trxframe(fio, fr);
 954             bFirst = FALSE;
 955             break;
 956         case efG96:
 957             /* Can not rewind a compressed file, so open it twice */
 958             if (!(*status)->persistent_line)
 959             {
 960                 /* allocate the persistent line */
 961                 snew((*status)->persistent_line, STRLEN+1);
 962             }
 963             read_g96_conf(gmx_fio_getfp(fio), fn, fr, (*status)->persistent_line);
 964             gmx_fio_close(fio);
 965             clear_trxframe(fr, FALSE);
 966             if (flags & (TRX_READ_X | TRX_NEED_X))
 967             {
 968                 snew(fr->x, fr->natoms);
 969             }
 970             if (flags & (TRX_READ_V | TRX_NEED_V))
 971             {
 972                 snew(fr->v, fr->natoms);
 973             }
 974             fio = (*status)->fio = gmx_fio_open(fn, "r");
 975             break;
 976         case efXTC:
 977             if (read_first_xtc(fio, &fr->natoms, &fr->step, &fr->time, fr->box, &fr->x,
 978                                &fr->prec, &bOK) == 0)
 979             {
 980                 assert(!bOK);
 981                 fr->not_ok = DATA_NOT_OK;
 982             }
 983             if (fr->not_ok)
 984             {
 985                 fr->natoms = 0;
 986                 printincomp(*status, fr);
 987             }
 988             else
 989             {
 990                 fr->bPrec = (fr->prec > 0);
 991                 fr->bStep = TRUE;
 992                 fr->bTime = TRUE;
 993                 fr->bX    = TRUE;
 994                 fr->bBox  = TRUE;
 995                 printcount(*status, oenv, fr->time, FALSE);
 996             }
 997             bFirst = FALSE;
 998             break;
 999         case efTNG:
1000             fr->step = -1;
1001             if (!gmx_read_next_tng_frame((*status)->tng, fr, NULL, 0))
1002             {
1003                 fr->not_ok = DATA_NOT_OK;
1004                 fr->natoms = 0;
1005                 printincomp(*status, fr);
1006             }
1007             else
1008             {
1009                 printcount(*status, oenv, fr->time, FALSE);
1010             }
1011             bFirst = FALSE;
1012             break;
1013         case efPDB:
1014             pdb_first_x(*status, gmx_fio_getfp(fio), fr);
1015             if (fr->natoms)
1016             {
1017                 printcount(*status, oenv, fr->time, FALSE);
1018             }
1019             bFirst = FALSE;
1020             break;
1021         case efGRO:
1022             if (gro_first_x_or_v(gmx_fio_getfp(fio), fr))
1023             {
1024                 printcount(*status, oenv, fr->time, FALSE);
1025             }
1026             bFirst = FALSE;
1027             break;
1028         default:
1029 #ifdef GMX_USE_PLUGINS
1030             fprintf(stderr, "The file format of %s is not a known trajectory format to GROMACS.\n"
1031                     "Please make sure that the file is a trajectory!\n"
1032                     "GROMACS will now assume it to be a trajectory and will try to open it using the VMD plug-ins.\n"
1033                     "This will only work in case the VMD plugins are found and it is a trajectory format supported by VMD.\n", fn);
1034             gmx_fio_fp_close(fio); /*only close the file without removing FIO entry*/
1035             if (!read_first_vmd_frame(fn, fr))
1036             {
1037                 gmx_fatal(FARGS, "Not supported in read_first_frame: %s", fn);
1038             }
1039 #else
1040             gmx_fatal(FARGS, "Not supported in read_first_frame: %s. Please make sure that the file is a trajectory.\n"
1041                       "GROMACS is not compiled with plug-in support. Thus it cannot read non-GROMACS trajectory formats using the VMD plug-ins.\n"
1042                       "Please compile with plug-in support if you want to read non-GROMACS trajectory formats.\n", fn);
1043 #endif
1044             break;
1045     }
1046
1047     /* Return FALSE if we read a frame that's past the set ending time. */
1048     if (!bFirst && (!(fr->flags & TRX_DONT_SKIP) && check_times(fr->time) > 0))
1049     {
1050         fr->t0 = fr->time;
1051         return FALSE;
1052     }
1053
1054     if (bFirst ||
1055         (!(fr->flags & TRX_DONT_SKIP) && check_times(fr->time) < 0))
1056     {
1057         /* Read a frame when no frame was read or the first was skipped */
1058         if (!read_next_frame(oenv, *status, fr))
1059         {
1060             return FALSE;
1061         }
1062     }
1063     fr->t0 = fr->time;
1064
1065     return (fr->natoms > 0);
1066 }
1067
1068 /***** C O O R D I N A T E   S T U F F *****/
1069
1070 int read_first_x(const output_env_t oenv, t_trxstatus **status, const char *fn,
1071                  real *t, rvec **x, matrix box)
1072 {
1073     t_trxframe fr;
1074
1075     read_first_frame(oenv, status, fn, &fr, TRX_NEED_X);
1076
1077     snew((*status)->xframe, 1);
1078     (*status)->nxframe   = 1;
1079     (*(*status)->xframe) = fr;
1080     *t                   = (*status)->xframe->time;
1081     *x                   = (*status)->xframe->x;
1082     copy_mat((*status)->xframe->box, box);
1083
1084     return (*status)->xframe->natoms;
1085 }
1086
1087 gmx_bool read_next_x(const output_env_t oenv, t_trxstatus *status, real *t,
1088                      rvec x[], matrix box)
1089 {
1090     gmx_bool bRet;
1091
1092     status->xframe->x = x;
1093     /*xframe[status].x = x;*/
1094     bRet = read_next_frame(oenv, status, status->xframe);
1095     *t   = status->xframe->time;
1096     copy_mat(status->xframe->box, box);
1097
1098     return bRet;
1099 }
1100
1101 void close_trj(t_trxstatus *status)
1102 {
1103     gmx_tng_close(&status->tng);
1104     if (status->fio)
1105     {
1106         gmx_fio_close(status->fio);
1107     }
1108
1109     /* The memory in status->xframe is lost here,
1110      * but the read_first_x/read_next_x functions are deprecated anyhow.
1111      * read_first_frame/read_next_frame and close_trx should be used.
1112      */
1113     sfree(status);
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, NULL, NULL, &natoms, NULL, NULL, NULL, top);
1132     if (ePBC)
1133     {
1134         *ePBC = epbc;
1135     }
1136
1137     return top;
1138 }