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 */
66 typedef enum { effXYZ, effXYZBox, effG87, effG87Box, effNR } eFileFormat;
78 char *persistent_line; /* Persistent line for reading g96 trajectories */
81 static void initcount(t_trxstatus *status)
86 static void status_init(t_trxstatus *status)
92 status->persistent_line=NULL;
96 int nframes_read(t_trxstatus *status)
98 return status->__frame;
101 static void printcount_(t_trxstatus *status, const output_env_t oenv,
102 const char *l,real t)
104 if ((status->__frame < 2*SKIP1 || status->__frame % SKIP1 == 0) &&
105 (status->__frame < 2*SKIP2 || status->__frame % SKIP2 == 0) &&
106 (status->__frame < 2*SKIP3 || status->__frame % SKIP3 == 0))
107 fprintf(stderr,"\r%-14s %6d time %8.3f ",l,status->__frame,
108 output_env_conv_time(oenv,t));
111 static void printcount(t_trxstatus *status, const output_env_t oenv,real t,
115 printcount_(status, oenv,bSkip ? "Skipping frame" : "Reading frame",t);
118 static void printlast(t_trxstatus *status, const output_env_t oenv,real t)
120 printcount_(status, oenv,"Last frame",t);
121 fprintf(stderr,"\n");
124 static void printincomp(t_trxstatus *status, t_trxframe *fr)
126 if (fr->not_ok & HEADER_NOT_OK)
127 fprintf(stderr,"WARNING: Incomplete header: nr %d time %g\n",
128 status->__frame+1,fr->time);
130 fprintf(stderr,"WARNING: Incomplete frame: nr %d time %g\n",
131 status->__frame+1,fr->time);
134 int prec2ndec(real prec)
137 gmx_fatal(FARGS,"DEATH HORROR prec (%g) <= 0 in prec2ndec",prec);
139 return (int)(log(prec)/log(10.0)+0.5);
143 t_fileio *trx_get_fileio(t_trxstatus *status)
150 void clear_trxframe(t_trxframe *fr,gmx_bool bFirst)
185 void set_trxframe_ePBC(t_trxframe *fr,int ePBC)
187 fr->bPBC = (ePBC == -1);
191 int write_trxframe_indexed(t_trxstatus *status,t_trxframe *fr,int nind,
192 atom_id *ind, gmx_conect gc)
195 rvec *xout=NULL,*vout=NULL,*fout=NULL;
204 switch (gmx_fio_getftp(status->fio)) {
210 gmx_fatal(FARGS,"Need coordinates to write a %s trajectory",
211 ftp2ext(gmx_fio_getftp(status->fio)));
215 switch (gmx_fio_getftp(status->fio)) {
220 for(i=0; i<nind; i++)
221 copy_rvec(fr->v[ind[i]],vout[i]);
225 for(i=0; i<nind; i++)
226 copy_rvec(fr->f[ind[i]],fout[i]);
233 for(i=0; i<nind; i++)
234 copy_rvec(fr->x[ind[i]],xout[i]);
241 switch (gmx_fio_getftp(status->fio)) {
243 write_xtc(status->fio,nind,fr->step,fr->time,fr->box,xout,prec);
247 fwrite_trn(status->fio,nframes_read(status),
248 fr->time,fr->step,fr->box,nind,xout,vout,fout);
255 gmx_fatal(FARGS,"Can not write a %s file without atom names",
256 ftp2ext(gmx_fio_getftp(status->fio)));
257 sprintf(title,"frame t= %.3f",fr->time);
258 if (gmx_fio_getftp(status->fio) == efGRO)
259 write_hconf_indexed_p(gmx_fio_getfp(status->fio),title,fr->atoms,nind,ind,
261 fr->x,fr->bV ? fr->v : NULL,fr->box);
263 write_pdbfile_indexed(gmx_fio_getfp(status->fio),title,fr->atoms,
264 fr->x,-1,fr->box,' ',fr->step,nind,ind,gc,TRUE);
267 write_gms(gmx_fio_getfp(status->fio),nind,xout,fr->box);
270 write_g96_conf(gmx_fio_getfp(status->fio),fr,nind,ind);
273 gmx_fatal(FARGS,"Sorry, write_trxframe_indexed can not write %s",
274 ftp2ext(gmx_fio_getftp(status->fio)));
278 switch (gmx_fio_getftp(status->fio)) {
282 if (vout) sfree(vout);
283 if (fout) sfree(fout);
296 int write_trxframe(t_trxstatus *status,t_trxframe *fr,gmx_conect gc)
306 switch (gmx_fio_getftp(status->fio)) {
312 gmx_fatal(FARGS,"Need coordinates to write a %s trajectory",
313 ftp2ext(gmx_fio_getftp(status->fio)));
317 switch (gmx_fio_getftp(status->fio)) {
319 write_xtc(status->fio,fr->natoms,fr->step,fr->time,fr->box,fr->x,prec);
323 fwrite_trn(status->fio,fr->step,fr->time,fr->lambda,fr->box,fr->natoms,
324 fr->bX ? fr->x:NULL,fr->bV ? fr->v:NULL ,fr->bF ? fr->f:NULL);
331 gmx_fatal(FARGS,"Can not write a %s file without atom names",
332 ftp2ext(gmx_fio_getftp(status->fio)));
333 sprintf(title,"frame t= %.3f",fr->time);
334 if (gmx_fio_getftp(status->fio) == efGRO)
335 write_hconf_p(gmx_fio_getfp(status->fio),title,fr->atoms,
336 prec2ndec(prec),fr->x,fr->bV ? fr->v : NULL,fr->box);
338 write_pdbfile(gmx_fio_getfp(status->fio),title,
339 fr->atoms,fr->x,fr->bPBC ? fr->ePBC : -1,fr->box,
340 ' ',fr->step,gc,TRUE);
343 write_gms(gmx_fio_getfp(status->fio),fr->natoms,fr->x,fr->box);
346 write_g96_conf(gmx_fio_getfp(status->fio),fr,-1,NULL);
349 gmx_fatal(FARGS,"Sorry, write_trxframe can not write %s",
350 ftp2ext(gmx_fio_getftp(status->fio)));
357 int write_trx(t_trxstatus *status,int nind,atom_id *ind,t_atoms *atoms,
358 int step,real time,matrix box,rvec x[],rvec *v,
363 clear_trxframe(&fr,TRUE);
368 fr.bAtoms = atoms!=NULL;
375 copy_mat(box,fr.box);
377 return write_trxframe_indexed(status,&fr,nind,ind,gc);
380 void close_trx(t_trxstatus *status)
382 gmx_fio_close(status->fio);
386 t_trxstatus *open_trx(const char *outfile,const char *filemode)
389 if (filemode[0]!='w' && filemode[0]!='a' && filemode[1]!='+')
390 gmx_fatal(FARGS,"Sorry, write_trx can only write");
395 stat->fio=gmx_fio_open(outfile,filemode);
399 static gmx_bool gmx_next_frame(t_trxstatus *status,t_trxframe *fr)
406 if (fread_trnheader(status->fio,&sh,&bOK)) {
407 fr->bDouble=sh.bDouble;
408 fr->natoms=sh.natoms;
414 fr->bFepState = TRUE;
415 fr->lambda = sh.lambda;
416 fr->bBox = sh.box_size>0;
417 if (fr->flags & (TRX_READ_X | TRX_NEED_X)) {
419 snew(fr->x,sh.natoms);
420 fr->bX = sh.x_size>0;
422 if (fr->flags & (TRX_READ_V | TRX_NEED_V)) {
424 snew(fr->v,sh.natoms);
425 fr->bV = sh.v_size>0;
427 if (fr->flags & (TRX_READ_F | TRX_NEED_F)) {
429 snew(fr->f,sh.natoms);
430 fr->bF = sh.f_size>0;
432 if (fread_htrn(status->fio,&sh,fr->box,fr->x,fr->v,fr->f))
435 fr->not_ok = DATA_NOT_OK;
438 fr->not_ok = HEADER_NOT_OK;
443 static void choose_file_format(FILE *fp)
451 printf(" Select File Format\n");
452 printf("---------------------------\n");
453 printf("1. XYZ File\n");
454 printf("2. XYZ File with Box\n");
455 printf("3. Gromos-87 Ascii Trajectory\n");
456 printf("4. Gromos-87 Ascii Trajectory with Box\n");
462 printf("\nChoice: ");
470 } while ((i < 0) || (i >= effNR));
473 stat->eFF = (eFileFormat) i;
475 for(m=0; (m<DIM); m++) stat->BOX[m]=0;
477 stat->bReadBox = (stat->eFF == effG87Box) || (stat->eFF == effXYZBox);
482 if( 5 != fscanf(fp,"%d%lf%lf%lf%lf",&stat->NATOMS,&stat->BOX[XX],&stat->BOX[YY],&stat->BOX[ZZ],&stat->DT))
484 gmx_fatal(FARGS,"Error reading natoms/box in file");
489 printf("GROMOS! OH DEAR...\n\n");
490 printf("Number of atoms ? ");
492 if (1 != scanf("%d",&stat->NATOMS))
494 gmx_fatal(FARGS,"Error reading natoms in file");
497 printf("Time between timeframes ? ");
499 if( 1 != scanf("%lf",&stat->DT))
501 gmx_fatal(FARGS,"Error reading dt from file");
504 if (stat->eFF == effG87) {
505 printf("Box X Y Z ? ");
507 if(3 != scanf("%lf%lf%lf",&stat->BOX[XX],&stat->BOX[YY],&stat->BOX[ZZ]))
509 gmx_fatal(FARGS,"Error reading box in file");
520 printf("Hellow World\n");
524 static gmx_bool do_read_xyz(t_trxstatus *status, FILE *fp,int natoms,
530 for(i=0; (i<natoms); i++) {
531 for(m=0; (m<DIM); m++) {
532 if (fscanf(fp,"%lf",&x0) != 1) {
534 fprintf(stderr,"error reading statusfile: x[%d][%d]\n",i,m);
541 if (status->bReadBox) {
542 for(m=0; (m<DIM); m++) {
543 if (fscanf(fp,"%lf",&x0) != 1)
551 static gmx_bool xyz_next_x(t_trxstatus *status, FILE *fp, const output_env_t oenv,
552 real *t, int natoms, rvec x[], matrix box)
553 /* Reads until a new x can be found (return TRUE)
554 * or eof (return FALSE)
560 while (!bTimeSet(TBEGIN) || (*t < rTimeValue(TBEGIN))) {
561 if (!do_read_xyz(status,fp,natoms,x,box))
563 printcount(status,oenv,*t,FALSE);
567 if (!bTimeSet(TEND) || (*t <= rTimeValue(TEND))) {
568 if (!do_read_xyz(status,fp,natoms,x,box)) {
569 printlast(status, oenv,*t);
572 printcount(status,oenv,*t,FALSE);
577 printlast(status,oenv,pt);
581 static int xyz_first_x(t_trxstatus *status, FILE *fp, const output_env_t oenv,
582 real *t, rvec **x, matrix box)
583 /* Reads fp, mallocs x, and returns x and box
584 * Returns natoms when successful, FALSE otherwise
592 choose_file_format(fp);
594 for(m=0; (m<DIM); m++)
595 box[m][m]=status->BOX[m];
597 snew(*x,status->NATOMS);
599 if (!xyz_next_x(status, fp,oenv,t,status->NATOMS,*x,box))
603 return status->NATOMS;
606 static gmx_bool pdb_next_x(t_trxstatus *status, FILE *fp,t_trxframe *fr)
610 int ePBC,model_nr,na;
611 char title[STRLEN],*time;
614 atoms.nr = fr->natoms;
617 /* the other pointers in atoms should not be accessed if these are NULL */
619 na=read_pdbfile(fp,title,&model_nr,&atoms,fr->x,&ePBC,boxpdb,TRUE,NULL);
620 set_trxframe_ePBC(fr,ePBC);
621 if (nframes_read(status)==0)
622 fprintf(stderr," '%s', %d atoms\n",title, fr->natoms);
626 fr->bBox = (boxpdb[XX][XX] != 0.0);
628 copy_mat(boxpdb,fr->box);
631 if (model_nr!=NOTSET) {
635 time=strstr(title," t= ");
638 sscanf(time+4,"%lf",&dbl);
642 /* this is a bit dirty, but it will work: if no time is read from
643 comment line in pdb file, set time to current frame number */
645 fr->time=(real)fr->step;
647 fr->time=(real)nframes_read(status);
652 if (na != fr->natoms)
653 gmx_fatal(FARGS,"Number of atoms in pdb frame %d is %d instead of %d",
654 nframes_read(status),na,fr->natoms);
659 static int pdb_first_x(t_trxstatus *status, FILE *fp, t_trxframe *fr)
663 fprintf(stderr,"Reading frames from pdb file");
665 get_pdb_coordnum(fp, &fr->natoms);
667 gmx_fatal(FARGS,"\nNo coordinates in pdb file\n");
669 snew(fr->x,fr->natoms);
670 pdb_next_x(status, fp, fr);
675 gmx_bool read_next_frame(const output_env_t oenv,t_trxstatus *status,t_trxframe *fr)
679 gmx_bool bOK,bRet,bMissingData=FALSE,bSkip=FALSE;
686 clear_trxframe(fr,FALSE);
690 switch (gmx_fio_getftp(status->fio)) {
693 bRet = gmx_next_frame(status,fr);
696 /* Checkpoint files can not contain mulitple frames */
699 read_g96_conf(gmx_fio_getfp(status->fio),NULL,fr,
700 status->persistent_line);
701 bRet = (fr->natoms > 0);
704 bRet = xyz_next_x(status, gmx_fio_getfp(status->fio),oenv,&fr->time,
705 fr->natoms, fr->x,fr->box);
712 * Sometimes is off by one frame
713 * and sometimes reports frame not present/file not seekable
715 /* DvdS 2005-05-31: this has been fixed along with the increased
716 * accuracy of the control over -b and -e options.
718 if (bTimeSet(TBEGIN) && (fr->time < rTimeValue(TBEGIN))) {
719 if (xtc_seek_time(status->fio, rTimeValue(TBEGIN),fr->natoms,TRUE)) {
720 gmx_fatal(FARGS,"Specified frame (time %f) doesn't exist or file corrupt/inconsistent.",
725 bRet = read_next_xtc(status->fio,fr->natoms,&fr->step,&fr->time,fr->box,
726 fr->x,&fr->prec,&bOK);
727 fr->bPrec = (bRet && fr->prec > 0);
733 /* Actually the header could also be not ok,
734 but from bOK from read_next_xtc this can't be distinguished */
735 fr->not_ok = DATA_NOT_OK;
739 bRet = pdb_next_x(status, gmx_fio_getfp(status->fio),fr);
742 bRet = gro_next_x_or_v(gmx_fio_getfp(status->fio),fr);
745 #ifdef GMX_USE_PLUGINS
746 bRet = read_next_vmd_frame(dummy,fr);
748 gmx_fatal(FARGS,"DEATH HORROR in read_next_frame ftp=%s,status=%s",
749 ftp2ext(gmx_fio_getftp(status->fio)),
750 gmx_fio_getname(status->fio));
755 bMissingData = (((fr->flags & TRX_NEED_X) && !fr->bX) ||
756 ((fr->flags & TRX_NEED_V) && !fr->bV) ||
757 ((fr->flags & TRX_NEED_F) && !fr->bF));
760 ct=check_times2(fr->time,fr->t0,fr->tpf,fr->tppf,fr->bDouble);
761 if (ct == 0 || ((fr->flags & TRX_DONT_SKIP) && ct<0)) {
762 printcount(status, oenv,fr->time,FALSE);
766 printcount(status, oenv,fr->time,TRUE);
772 } while (bRet && (bMissingData || bSkip));
775 printlast(status, oenv,pt);
777 printincomp(status, fr);
783 int read_first_frame(const output_env_t oenv,t_trxstatus **status,
784 const char *fn,t_trxframe *fr,int flags)
790 clear_trxframe(fr,TRUE);
797 status_init( *status );
798 (*status)->nxframe=1;
801 fio = (*status)->fio =gmx_fio_open(fn,"r");
802 switch (gmx_fio_getftp(fio))
808 read_checkpoint_trxframe(fio,fr);
812 /* Can not rewind a compressed file, so open it twice */
813 if (!(*status)->persistent_line)
815 /* allocate the persistent line */
816 snew((*status)->persistent_line, STRLEN+1);
818 read_g96_conf(gmx_fio_getfp(fio),fn,fr, (*status)->persistent_line);
820 clear_trxframe(fr,FALSE);
821 if (flags & (TRX_READ_X | TRX_NEED_X))
822 snew(fr->x,fr->natoms);
823 if (flags & (TRX_READ_V | TRX_NEED_V))
824 snew(fr->v,fr->natoms);
825 fio = (*status)->fio =gmx_fio_open(fn,"r");
828 fr->natoms=xyz_first_x(*status, gmx_fio_getfp(fio),oenv,&fr->time,
834 printcount(*status,oenv,fr->time,FALSE);
839 if (read_first_xtc(fio,&fr->natoms,&fr->step,&fr->time,fr->box,&fr->x,
840 &fr->prec,&bOK) == 0) {
842 gmx_fatal(FARGS,"No XTC!\n");
844 fr->not_ok = DATA_NOT_OK;
849 printincomp(*status,fr);
851 fr->bPrec = (fr->prec > 0);
856 printcount(*status,oenv,fr->time,FALSE);
861 pdb_first_x(*status, gmx_fio_getfp(fio),fr);
863 printcount(*status,oenv,fr->time,FALSE);
867 if (gro_first_x_or_v(gmx_fio_getfp(fio),fr))
868 printcount(*status,oenv,fr->time,FALSE);
872 #ifdef GMX_USE_PLUGINS
873 fprintf(stderr,"The file format of %s is not a known trajectory format to GROMACS.\n"
874 "Please make sure that the file is a trajectory!\n"
875 "GROMACS will now assume it to be a trajectory and will try to open it using the VMD plug-ins.\n"
876 "This will only work in case the VMD plugins are found and it is a trajectory format supported by VMD.\n",fn);
877 gmx_fio_fp_close(fio); /*only close the file without removing FIO entry*/
878 if (!read_first_vmd_frame(&dummy,fn,fr,flags))
880 gmx_fatal(FARGS,"Not supported in read_first_frame: %s",fn);
883 gmx_fatal(FARGS,"Not supported in read_first_frame: %s. Please make sure that the file is a trajectory.\n"
884 "GROMACS is not compiled with plug-in support. Thus it cannot read non-GROMACS trajectory formats using the VMD plug-ins.\n"
885 "Please compile with plug-in support if you want to read non-GROMACS trajectory formats.\n",fn);
890 /* Return FALSE if we read a frame that's past the set ending time. */
891 if (!bFirst && (!(fr->flags & TRX_DONT_SKIP) && check_times(fr->time) > 0)) {
897 (!(fr->flags & TRX_DONT_SKIP) && check_times(fr->time) < 0))
898 /* Read a frame when no frame was read or the first was skipped */
899 if (!read_next_frame(oenv,*status,fr))
903 return (fr->natoms > 0);
906 /***** C O O R D I N A T E S T U F F *****/
908 int read_first_x(const output_env_t oenv,t_trxstatus **status,const char *fn,
909 real *t,rvec **x,matrix box)
913 read_first_frame(oenv,status,fn,&fr,TRX_NEED_X);
915 snew((*status)->xframe, 1);
916 (*status)->nxframe=1;
917 (*(*status)->xframe) = fr;
918 *t = (*status)->xframe->time;
919 *x = (*status)->xframe->x;
920 copy_mat((*status)->xframe->box,box);
922 return (*status)->xframe->natoms;
925 gmx_bool read_next_x(const output_env_t oenv, t_trxstatus *status,real *t,
926 int natoms, rvec x[], matrix box)
930 status->xframe->x= x;
931 /*xframe[status].x = x;*/
932 bRet = read_next_frame(oenv,status,status->xframe);
933 *t = status->xframe->time;
934 copy_mat(status->xframe->box,box);
939 void close_trj(t_trxstatus *status)
941 gmx_fio_close(status->fio);
942 /* The memory in status->xframe is lost here,
943 * but the read_first_x/read_next_x functions are deprecated anyhow.
944 * read_first_frame/read_next_frame and close_trx should be used.
949 void rewind_trj(t_trxstatus *status)
953 gmx_fio_rewind(status->fio);
956 /***** V E L O C I T Y S T U F F *****/
958 static void clear_v(t_trxframe *fr)
963 for(i=0; i<fr->natoms; i++)
964 clear_rvec(fr->v[i]);