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