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