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