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