3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code forxd
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * GROningen Mixture of Alchemy and Childrens' Stories
56 #include "checkpoint.h"
60 /* defines for frame counter output */
65 /* Globals for gromos-87 input */
67 effXYZ, effXYZBox, effG87, effG87Box, effNR
80 char *persistent_line; /* Persistent line for reading g96 trajectories */
83 static void initcount(t_trxstatus *status)
88 static void status_init(t_trxstatus *status)
91 status->xframe = NULL;
94 status->persistent_line = NULL;
98 int nframes_read(t_trxstatus *status)
100 return status->__frame;
103 static void printcount_(t_trxstatus *status, const output_env_t oenv,
104 const char *l, real t)
106 if ((status->__frame < 2*SKIP1 || status->__frame % SKIP1 == 0) &&
107 (status->__frame < 2*SKIP2 || status->__frame % SKIP2 == 0) &&
108 (status->__frame < 2*SKIP3 || status->__frame % SKIP3 == 0))
110 fprintf(stderr, "\r%-14s %6d time %8.3f ", l, status->__frame,
111 output_env_conv_time(oenv, t));
115 static void printcount(t_trxstatus *status, const output_env_t oenv, real t,
119 printcount_(status, oenv, bSkip ? "Skipping frame" : "Reading frame", t);
122 static void printlast(t_trxstatus *status, const output_env_t oenv, real t)
124 printcount_(status, oenv, "Last frame", t);
125 fprintf(stderr, "\n");
128 static void printincomp(t_trxstatus *status, t_trxframe *fr)
130 if (fr->not_ok & HEADER_NOT_OK)
132 fprintf(stderr, "WARNING: Incomplete header: nr %d time %g\n",
133 status->__frame+1, fr->time);
137 fprintf(stderr, "WARNING: Incomplete frame: nr %d time %g\n",
138 status->__frame+1, fr->time);
142 int prec2ndec(real prec)
146 gmx_fatal(FARGS, "DEATH HORROR prec (%g) <= 0 in prec2ndec", prec);
149 return (int)(log(prec)/log(10.0)+0.5);
153 t_fileio *trx_get_fileio(t_trxstatus *status)
160 void clear_trxframe(t_trxframe *fr, gmx_bool bFirst)
196 void set_trxframe_ePBC(t_trxframe *fr, int ePBC)
198 fr->bPBC = (ePBC == -1);
202 int write_trxframe_indexed(t_trxstatus *status, t_trxframe *fr, int nind,
203 const atom_id *ind, gmx_conect gc)
206 rvec *xout = NULL, *vout = NULL, *fout = NULL;
219 switch (gmx_fio_getftp(status->fio))
227 gmx_fatal(FARGS, "Need coordinates to write a %s trajectory",
228 ftp2ext(gmx_fio_getftp(status->fio)));
233 switch (gmx_fio_getftp(status->fio))
240 for (i = 0; i < nind; i++)
242 copy_rvec(fr->v[ind[i]], vout[i]);
248 for (i = 0; i < nind; i++)
250 copy_rvec(fr->f[ind[i]], fout[i]);
259 for (i = 0; i < nind; i++)
261 copy_rvec(fr->x[ind[i]], xout[i]);
269 switch (gmx_fio_getftp(status->fio))
272 write_xtc(status->fio, nind, fr->step, fr->time, fr->box, xout, prec);
276 fwrite_trn(status->fio, nframes_read(status),
277 fr->time, fr->step, fr->box, nind, xout, vout, fout);
285 gmx_fatal(FARGS, "Can not write a %s file without atom names",
286 ftp2ext(gmx_fio_getftp(status->fio)));
288 sprintf(title, "frame t= %.3f", fr->time);
289 if (gmx_fio_getftp(status->fio) == efGRO)
291 write_hconf_indexed_p(gmx_fio_getfp(status->fio), title, fr->atoms, nind, ind,
293 fr->x, fr->bV ? fr->v : NULL, fr->box);
297 write_pdbfile_indexed(gmx_fio_getfp(status->fio), title, fr->atoms,
298 fr->x, -1, fr->box, ' ', fr->step, nind, ind, gc, TRUE);
302 write_gms(gmx_fio_getfp(status->fio), nind, xout, fr->box);
305 write_g96_conf(gmx_fio_getfp(status->fio), fr, nind, ind);
308 gmx_fatal(FARGS, "Sorry, write_trxframe_indexed can not write %s",
309 ftp2ext(gmx_fio_getftp(status->fio)));
313 switch (gmx_fio_getftp(status->fio))
338 int write_trxframe(t_trxstatus *status, t_trxframe *fr, gmx_conect gc)
352 switch (gmx_fio_getftp(status->fio))
360 gmx_fatal(FARGS, "Need coordinates to write a %s trajectory",
361 ftp2ext(gmx_fio_getftp(status->fio)));
366 switch (gmx_fio_getftp(status->fio))
369 write_xtc(status->fio, fr->natoms, fr->step, fr->time, fr->box, fr->x, prec);
373 fwrite_trn(status->fio, fr->step, fr->time, fr->lambda, fr->box, fr->natoms,
374 fr->bX ? fr->x : NULL, fr->bV ? fr->v : NULL, fr->bF ? fr->f : NULL);
382 gmx_fatal(FARGS, "Can not write a %s file without atom names",
383 ftp2ext(gmx_fio_getftp(status->fio)));
385 sprintf(title, "frame t= %.3f", fr->time);
386 if (gmx_fio_getftp(status->fio) == efGRO)
388 write_hconf_p(gmx_fio_getfp(status->fio), title, fr->atoms,
389 prec2ndec(prec), fr->x, fr->bV ? fr->v : NULL, fr->box);
393 write_pdbfile(gmx_fio_getfp(status->fio), title,
394 fr->atoms, fr->x, fr->bPBC ? fr->ePBC : -1, fr->box,
395 ' ', fr->step, gc, TRUE);
399 write_gms(gmx_fio_getfp(status->fio), fr->natoms, fr->x, fr->box);
402 write_g96_conf(gmx_fio_getfp(status->fio), fr, -1, NULL);
405 gmx_fatal(FARGS, "Sorry, write_trxframe can not write %s",
406 ftp2ext(gmx_fio_getftp(status->fio)));
413 int write_trx(t_trxstatus *status, int nind, const atom_id *ind, t_atoms *atoms,
414 int step, real time, matrix box, rvec x[], rvec *v,
419 clear_trxframe(&fr, TRUE);
424 fr.bAtoms = atoms != NULL;
431 copy_mat(box, fr.box);
433 return write_trxframe_indexed(status, &fr, nind, ind, gc);
436 void close_trx(t_trxstatus *status)
438 gmx_fio_close(status->fio);
442 t_trxstatus *open_trx(const char *outfile, const char *filemode)
445 if (filemode[0] != 'w' && filemode[0] != 'a' && filemode[1] != '+')
447 gmx_fatal(FARGS, "Sorry, write_trx can only write");
453 stat->fio = gmx_fio_open(outfile, filemode);
457 static gmx_bool gmx_next_frame(t_trxstatus *status, t_trxframe *fr)
464 if (fread_trnheader(status->fio, &sh, &bOK))
466 fr->bDouble = sh.bDouble;
467 fr->natoms = sh.natoms;
473 fr->bFepState = TRUE;
474 fr->lambda = sh.lambda;
475 fr->bBox = sh.box_size > 0;
476 if (fr->flags & (TRX_READ_X | TRX_NEED_X))
480 snew(fr->x, sh.natoms);
482 fr->bX = sh.x_size > 0;
484 if (fr->flags & (TRX_READ_V | TRX_NEED_V))
488 snew(fr->v, sh.natoms);
490 fr->bV = sh.v_size > 0;
492 if (fr->flags & (TRX_READ_F | TRX_NEED_F))
496 snew(fr->f, sh.natoms);
498 fr->bF = sh.f_size > 0;
500 if (fread_htrn(status->fio, &sh, fr->box, fr->x, fr->v, fr->f))
506 fr->not_ok = DATA_NOT_OK;
512 fr->not_ok = HEADER_NOT_OK;
518 static void choose_file_format(FILE *fp)
526 printf(" Select File Format\n");
527 printf("---------------------------\n");
528 printf("1. XYZ File\n");
529 printf("2. XYZ File with Box\n");
530 printf("3. Gromos-87 Ascii Trajectory\n");
531 printf("4. Gromos-87 Ascii Trajectory with Box\n");
538 printf("\nChoice: ");
542 rc = scanf("%d", &i);
547 while ((i < 0) || (i >= effNR));
550 stat->eFF = (eFileFormat) i;
552 for (m = 0; (m < DIM); m++)
557 stat->bReadBox = (stat->eFF == effG87Box) || (stat->eFF == effXYZBox);
563 if (5 != fscanf(fp, "%d%lf%lf%lf%lf", &stat->NATOMS, &stat->BOX[XX], &stat->BOX[YY], &stat->BOX[ZZ], &stat->DT))
565 gmx_fatal(FARGS, "Error reading natoms/box in file");
570 printf("GROMOS! OH DEAR...\n\n");
571 printf("Number of atoms ? ");
573 if (1 != scanf("%d", &stat->NATOMS))
575 gmx_fatal(FARGS, "Error reading natoms in file");
578 printf("Time between timeframes ? ");
580 if (1 != scanf("%lf", &stat->DT))
582 gmx_fatal(FARGS, "Error reading dt from file");
585 if (stat->eFF == effG87)
587 printf("Box X Y Z ? ");
589 if (3 != scanf("%lf%lf%lf", &stat->BOX[XX], &stat->BOX[YY], &stat->BOX[ZZ]))
591 gmx_fatal(FARGS, "Error reading box in file");
604 printf("Hellow World\n");
608 static gmx_bool do_read_xyz(t_trxstatus *status, FILE *fp, int natoms,
609 rvec x[], matrix box)
614 for (i = 0; (i < natoms); i++)
616 for (m = 0; (m < DIM); m++)
618 if (fscanf(fp, "%lf", &x0) != 1)
622 fprintf(stderr, "error reading statusfile: x[%d][%d]\n", i, m);
630 if (status->bReadBox)
632 for (m = 0; (m < DIM); m++)
634 if (fscanf(fp, "%lf", &x0) != 1)
644 static gmx_bool xyz_next_x(t_trxstatus *status, FILE *fp, const output_env_t oenv,
645 real *t, int natoms, rvec x[], matrix box)
646 /* Reads until a new x can be found (return TRUE)
647 * or eof (return FALSE)
653 while (!bTimeSet(TBEGIN) || (*t < rTimeValue(TBEGIN)))
655 if (!do_read_xyz(status, fp, natoms, x, box))
659 printcount(status, oenv, *t, FALSE);
663 if (!bTimeSet(TEND) || (*t <= rTimeValue(TEND)))
665 if (!do_read_xyz(status, fp, natoms, x, box))
667 printlast(status, oenv, *t);
670 printcount(status, oenv, *t, FALSE);
675 printlast(status, oenv, pt);
679 static int xyz_first_x(t_trxstatus *status, FILE *fp, const output_env_t oenv,
680 real *t, rvec **x, matrix box)
681 /* Reads fp, mallocs x, and returns x and box
682 * Returns natoms when successful, FALSE otherwise
690 choose_file_format(fp);
692 for (m = 0; (m < DIM); m++)
694 box[m][m] = status->BOX[m];
697 snew(*x, status->NATOMS);
699 if (!xyz_next_x(status, fp, oenv, t, status->NATOMS, *x, box))
705 return status->NATOMS;
708 static gmx_bool pdb_next_x(t_trxstatus *status, FILE *fp, t_trxframe *fr)
712 int ePBC, model_nr, na;
713 char title[STRLEN], *time;
716 atoms.nr = fr->natoms;
718 atoms.pdbinfo = NULL;
719 /* the other pointers in atoms should not be accessed if these are NULL */
721 na = read_pdbfile(fp, title, &model_nr, &atoms, fr->x, &ePBC, boxpdb, TRUE, NULL);
722 set_trxframe_ePBC(fr, ePBC);
723 if (nframes_read(status) == 0)
725 fprintf(stderr, " '%s', %d atoms\n", title, fr->natoms);
730 fr->bBox = (boxpdb[XX][XX] != 0.0);
733 copy_mat(boxpdb, fr->box);
736 if (model_nr != NOTSET)
741 time = strstr(title, " t= ");
745 sscanf(time+4, "%lf", &dbl);
746 fr->time = (real)dbl;
751 /* this is a bit dirty, but it will work: if no time is read from
752 comment line in pdb file, set time to current frame number */
755 fr->time = (real)fr->step;
759 fr->time = (real)nframes_read(status);
768 if (na != fr->natoms)
770 gmx_fatal(FARGS, "Number of atoms in pdb frame %d is %d instead of %d",
771 nframes_read(status), na, fr->natoms);
777 static int pdb_first_x(t_trxstatus *status, FILE *fp, t_trxframe *fr)
781 fprintf(stderr, "Reading frames from pdb file");
783 get_pdb_coordnum(fp, &fr->natoms);
786 gmx_fatal(FARGS, "\nNo coordinates in pdb file\n");
789 snew(fr->x, fr->natoms);
790 pdb_next_x(status, fp, fr);
795 gmx_bool read_next_frame(const output_env_t oenv, t_trxstatus *status, t_trxframe *fr)
799 gmx_bool bOK, bRet, bMissingData = FALSE, bSkip = FALSE;
807 clear_trxframe(fr, FALSE);
811 switch (gmx_fio_getftp(status->fio))
815 bRet = gmx_next_frame(status, fr);
818 /* Checkpoint files can not contain mulitple frames */
821 read_g96_conf(gmx_fio_getfp(status->fio), NULL, fr,
822 status->persistent_line);
823 bRet = (fr->natoms > 0);
826 bRet = xyz_next_x(status, gmx_fio_getfp(status->fio), oenv, &fr->time,
827 fr->natoms, fr->x, fr->box);
834 * Sometimes is off by one frame
835 * and sometimes reports frame not present/file not seekable
837 /* DvdS 2005-05-31: this has been fixed along with the increased
838 * accuracy of the control over -b and -e options.
840 if (bTimeSet(TBEGIN) && (fr->time < rTimeValue(TBEGIN)))
842 if (xtc_seek_time(status->fio, rTimeValue(TBEGIN), fr->natoms, TRUE))
844 gmx_fatal(FARGS, "Specified frame (time %f) doesn't exist or file corrupt/inconsistent.",
849 bRet = read_next_xtc(status->fio, fr->natoms, &fr->step, &fr->time, fr->box,
850 fr->x, &fr->prec, &bOK);
851 fr->bPrec = (bRet && fr->prec > 0);
858 /* Actually the header could also be not ok,
859 but from bOK from read_next_xtc this can't be distinguished */
860 fr->not_ok = DATA_NOT_OK;
864 bRet = pdb_next_x(status, gmx_fio_getfp(status->fio), fr);
867 bRet = gro_next_x_or_v(gmx_fio_getfp(status->fio), fr);
870 #ifdef GMX_USE_PLUGINS
871 bRet = read_next_vmd_frame(dummy, fr);
873 gmx_fatal(FARGS, "DEATH HORROR in read_next_frame ftp=%s,status=%s",
874 ftp2ext(gmx_fio_getftp(status->fio)),
875 gmx_fio_getname(status->fio));
881 bMissingData = (((fr->flags & TRX_NEED_X) && !fr->bX) ||
882 ((fr->flags & TRX_NEED_V) && !fr->bV) ||
883 ((fr->flags & TRX_NEED_F) && !fr->bF));
887 ct = check_times2(fr->time, fr->t0, fr->tpf, fr->tppf, fr->bDouble);
888 if (ct == 0 || ((fr->flags & TRX_DONT_SKIP) && ct < 0))
890 printcount(status, oenv, fr->time, FALSE);
898 printcount(status, oenv, fr->time, TRUE);
905 while (bRet && (bMissingData || bSkip));
909 printlast(status, oenv, pt);
912 printincomp(status, fr);
919 int read_first_frame(const output_env_t oenv, t_trxstatus **status,
920 const char *fn, t_trxframe *fr, int flags)
923 gmx_bool bFirst, bOK;
926 clear_trxframe(fr, TRUE);
933 status_init( *status );
934 (*status)->nxframe = 1;
937 fio = (*status)->fio = gmx_fio_open(fn, "r");
938 switch (gmx_fio_getftp(fio))
944 read_checkpoint_trxframe(fio, fr);
948 /* Can not rewind a compressed file, so open it twice */
949 if (!(*status)->persistent_line)
951 /* allocate the persistent line */
952 snew((*status)->persistent_line, STRLEN+1);
954 read_g96_conf(gmx_fio_getfp(fio), fn, fr, (*status)->persistent_line);
956 clear_trxframe(fr, FALSE);
957 if (flags & (TRX_READ_X | TRX_NEED_X))
959 snew(fr->x, fr->natoms);
961 if (flags & (TRX_READ_V | TRX_NEED_V))
963 snew(fr->v, fr->natoms);
965 fio = (*status)->fio = gmx_fio_open(fn, "r");
968 fr->natoms = xyz_first_x(*status, gmx_fio_getfp(fio), oenv, &fr->time,
975 printcount(*status, oenv, fr->time, FALSE);
980 if (read_first_xtc(fio, &fr->natoms, &fr->step, &fr->time, fr->box, &fr->x,
981 &fr->prec, &bOK) == 0)
985 gmx_fatal(FARGS, "No XTC!\n");
989 fr->not_ok = DATA_NOT_OK;
995 printincomp(*status, fr);
999 fr->bPrec = (fr->prec > 0);
1004 printcount(*status, oenv, fr->time, FALSE);
1009 pdb_first_x(*status, gmx_fio_getfp(fio), fr);
1012 printcount(*status, oenv, fr->time, FALSE);
1017 if (gro_first_x_or_v(gmx_fio_getfp(fio), fr))
1019 printcount(*status, oenv, fr->time, FALSE);
1024 #ifdef GMX_USE_PLUGINS
1025 fprintf(stderr, "The file format of %s is not a known trajectory format to GROMACS.\n"
1026 "Please make sure that the file is a trajectory!\n"
1027 "GROMACS will now assume it to be a trajectory and will try to open it using the VMD plug-ins.\n"
1028 "This will only work in case the VMD plugins are found and it is a trajectory format supported by VMD.\n", fn);
1029 gmx_fio_fp_close(fio); /*only close the file without removing FIO entry*/
1030 if (!read_first_vmd_frame(&dummy, fn, fr, flags))
1032 gmx_fatal(FARGS, "Not supported in read_first_frame: %s", fn);
1035 gmx_fatal(FARGS, "Not supported in read_first_frame: %s. Please make sure that the file is a trajectory.\n"
1036 "GROMACS is not compiled with plug-in support. Thus it cannot read non-GROMACS trajectory formats using the VMD plug-ins.\n"
1037 "Please compile with plug-in support if you want to read non-GROMACS trajectory formats.\n", fn);
1042 /* Return FALSE if we read a frame that's past the set ending time. */
1043 if (!bFirst && (!(fr->flags & TRX_DONT_SKIP) && check_times(fr->time) > 0))
1050 (!(fr->flags & TRX_DONT_SKIP) && check_times(fr->time) < 0))
1052 /* Read a frame when no frame was read or the first was skipped */
1053 if (!read_next_frame(oenv, *status, fr))
1060 return (fr->natoms > 0);
1063 /***** C O O R D I N A T E S T U F F *****/
1065 int read_first_x(const output_env_t oenv, t_trxstatus **status, const char *fn,
1066 real *t, rvec **x, matrix box)
1070 read_first_frame(oenv, status, fn, &fr, TRX_NEED_X);
1072 snew((*status)->xframe, 1);
1073 (*status)->nxframe = 1;
1074 (*(*status)->xframe) = fr;
1075 *t = (*status)->xframe->time;
1076 *x = (*status)->xframe->x;
1077 copy_mat((*status)->xframe->box, box);
1079 return (*status)->xframe->natoms;
1082 gmx_bool read_next_x(const output_env_t oenv, t_trxstatus *status, real *t,
1083 int natoms, rvec x[], matrix box)
1087 status->xframe->x = x;
1088 /*xframe[status].x = x;*/
1089 bRet = read_next_frame(oenv, status, status->xframe);
1090 *t = status->xframe->time;
1091 copy_mat(status->xframe->box, box);
1096 void close_trj(t_trxstatus *status)
1098 gmx_fio_close(status->fio);
1099 /* The memory in status->xframe is lost here,
1100 * but the read_first_x/read_next_x functions are deprecated anyhow.
1101 * read_first_frame/read_next_frame and close_trx should be used.
1106 void rewind_trj(t_trxstatus *status)
1110 gmx_fio_rewind(status->fio);