src/gromacs/fileio/trxio.c

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