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