src/gmxlib/enxio.c

   1 /*
   2  * $Id$
   3  *
   4  *       This source code is part of
   5  *
   6  *        G   R   O   M   A   C   S
   7  *
   8  * GROningen MAchine for Chemical Simulations
   9  *
  10  *               VERSION 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * Green Red Orange Magenta Azure Cyan Skyblue
  28  */
  29 static char *SRCID_enxio_c = "$Id$";
  30
  31 #include "futil.h"
  32 #include "string2.h"
  33 #include "fatal.h"
  34 #include "smalloc.h"
  35 #include "gmxfio.h"
  36 #include "enxio.h"
  37
  38 typedef struct {
  39   real t;               /* Time frame                           */
  40   int  step;            /* MD step                              */
  41   int  nre;             /* Number of energies                   */
  42   int  ndisre;          /* Number of disre blocks               */
  43   int  nuser;           /* User definable number                */
  44   int  e_size;          /* Size (in bytes) of energies          */
  45   int  d_size;          /* Size (in bytes) of disre blocks      */
  46   int  u_size;          /* Size (in bytes) of user blocks       */
  47 } t_eheader;
  48
  49 static void wr_ener_nms(FILE *out,int nre,char *nms[])
  50 {
  51   int i;
  52
  53   fprintf(out,"%d\n",nre);
  54   for(i=0; (i<nre); i++)
  55     fprintf(out,"%s\n",nms[i]);
  56 }
  57
  58 static void rd_ener_nms(FILE *in,int *nre,char ***nm)
  59 {
  60   char line[256];
  61   int  i;
  62
  63   fgets2(line,255,in);
  64   if (sscanf(line,"%d",nre) == 0) {
  65     *nre=0;
  66     return;
  67   }
  68   snew((*nm),*nre);
  69   for(i=0; (i< (*nre)); i++) {
  70     fgets2(line,255,in);
  71     trim(line);
  72     (*nm)[i]=strdup(line);
  73   }
  74 }
  75
  76 static void edr_nms(XDR *xdr,int *nre,char ***nms)
  77 {
  78   int i;
  79   char **NM;
  80
  81   NM=*nms;
  82
  83   if (!xdr_int(xdr,nre)) {
  84     *nre=0;
  85     return;
  86   }
  87   if (NM == NULL) {
  88     snew(NM,*nre);
  89   }
  90   for(i=0; (i<*nre); i++)
  91     xdr_string(xdr,&(NM[i]),NM[i] ? strlen(NM[i]) : STRLEN);
  92   *nms=NM;
  93 }
  94
  95 static bool do_eheader(int fp,t_eheader *eh,bool *bOK)
  96 {
  97   bool bRead = fio_getread(fp);
  98
  99   *bOK=TRUE;
 100   if (!do_real(eh->t))      return FALSE;
 101   if (!do_int (eh->step))   *bOK = FALSE;
 102   if (!do_int (eh->nre))    *bOK = FALSE;
 103   if (!do_int (eh->ndisre)) *bOK = FALSE;
 104   if (!do_int (eh->nuser))  *bOK = FALSE;
 105   if (!do_int (eh->e_size)) *bOK = FALSE;
 106   if (!do_int (eh->d_size)) *bOK = FALSE;
 107   if (!do_int (eh->u_size)) *bOK = FALSE;
 108
 109   return *bOK;
 110 }
 111
 112 void do_enxnms(int fp,int *nre,char ***nms)
 113 {
 114   bool bRead;
 115
 116   bRead = fio_getread(fp);
 117   if (fio_getftp(fp) == efEDR) {
 118     fio_select(fp);
 119     edr_nms(fio_getxdr(fp),nre,nms);
 120   }
 121   else if (bRead)
 122     rd_ener_nms(fio_getfp(fp),nre,nms);
 123   else
 124     wr_ener_nms(fio_getfp(fp),*nre,*nms);
 125 }
 126
 127 void close_enx(int fp)
 128 {
 129   fio_close(fp);
 130 }
 131
 132 static bool empty_file(char *fn)
 133 {
 134   FILE *fp;
 135   char dum;
 136   bool bEmpty;
 137
 138   fp = ffopen(fn,"r");
 139   fread(&dum,sizeof(dum),1,fp);
 140   bEmpty = feof(fp);
 141   fclose(fp);
 142
 143   return bEmpty;
 144 }
 145
 146 static int framenr;
 147
 148 int open_enx(char *fn,char *mode)
 149 {
 150   int       i,fp;
 151   char      **nm=NULL;
 152   int       nre;
 153   t_eheader *eh;
 154   bool bDum=TRUE;
 155
 156   if (mode[0]=='r') {
 157     fp=fio_open(fn,mode);
 158     fio_select(fp);
 159     fio_setprecision(fp,FALSE);
 160     do_enxnms(fp,&nre,&nm);
 161     snew(eh,1);
 162     do_eheader(fp,eh,&bDum);
 163
 164     /* Now check whether this file is in single precision */
 165     if (((eh->e_size && (eh->nre == nre) && (nre*4*sizeof(float) == eh->e_size)) ||
 166          (eh->d_size && (eh->ndisre*sizeof(float)*2+sizeof(int) == eh->d_size)))){
 167       fprintf(stderr,"Opened %s as single precision energy file\n",fn);
 168       for(i=0; (i<nre); i++)
 169         sfree(nm[i]);
 170       sfree(nm);
 171     }
 172     else {
 173       fio_rewind(fp);
 174       fio_select(fp);
 175       fio_setprecision(fp,TRUE);
 176       do_enxnms(fp,&nre,&nm);
 177       do_eheader(fp,eh,&bDum);
 178       if (((eh->e_size && (eh->nre == nre) && (nre*4*sizeof(double) == eh->e_size)) ||
 179            (eh->d_size && (eh->ndisre*sizeof(double)*2+sizeof(int) == eh->d_size))))
 180         fprintf(stderr,"Opened %s as double precision energy file\n",fn);
 181       else {
 182         if (empty_file(fn))
 183           fatal_error(0,"File %s is empty",fn);
 184         else
 185           fatal_error(0,"Energy file %s not recognized, maybe different CPU?",
 186                       fn);
 187       }
 188       for(i=0; (i<nre); i++)
 189           sfree(nm[i]);
 190       sfree(nm);
 191     }
 192     sfree(eh);
 193     fio_rewind(fp);
 194   }
 195   else
 196     fp = fio_open(fn,mode);
 197
 198   framenr=0;
 199
 200   return fp;
 201 }
 202
 203 bool do_enx(int fp,real *t,int *step,int *nre,t_energy ener[],
 204             int *ndr,t_drblock *drblock)
 205 {
 206   int       i;
 207   t_eheader eh;
 208   bool      bRead,bOK,bOK1;
 209   real      tmp;
 210
 211   bOK = TRUE;
 212   bRead = fio_getread(fp);
 213   if (!bRead) {
 214     eh.t      = *t;
 215     eh.step   = *step;
 216     eh.nre    = *nre;
 217     eh.ndisre = drblock ? (drblock->ndr) : 0;
 218     eh.nuser  = 0;
 219     eh.e_size = (*nre)*sizeof(ener[0].e)*4;
 220     eh.d_size = eh.ndisre ?
 221       (drblock->ndr*(sizeof(drblock->rav[0])+sizeof(drblock->rt[0]))) : 0;
 222     eh.u_size = 0;
 223   }
 224   fio_select(fp);
 225
 226   if (!do_eheader(fp,&eh,&bOK)) {
 227     if (bRead) {
 228         fprintf(stderr,"\rLast frame read %d                          ",framenr-1);
 229         if (!bOK)
 230           fprintf(stderr,"\nWARNING: Incomplete frame: nr %6d time %8.3f\n",
 231                   framenr,eh.t);
 232     }
 233     return FALSE;
 234   }
 235   if (bRead) {
 236     if ( ( framenr<10 ) || ( framenr%10 == 0) )
 237       fprintf(stderr,"\rReading frame %6d time %8.3f           ",framenr,eh.t);
 238     framenr++;
 239   }
 240   *t    = eh.t;
 241   *step = eh.step;
 242   *nre  = eh.nre;
 243   for(i=0; (i<*nre); i++) {
 244     bOK = bOK && do_real(ener[i].e);
 245     tmp=ener[i].eav;
 246     if((tmp/(*step+1))<GMX_REAL_EPS)
 247       tmp=0;
 248     bOK = bOK && do_real(tmp);
 249     ener[i].eav=tmp;
 250     tmp=ener[i].esum;
 251     bOK = bOK && do_real(tmp);
 252     ener[i].esum=tmp;
 253     bOK = bOK && do_real(ener[i].e2sum);
 254   }
 255   if (eh.ndisre) {
 256     if (drblock->ndr == 0) {
 257       snew(drblock->rav,eh.ndisre);
 258       snew(drblock->rt,eh.ndisre);
 259       drblock->ndr = eh.ndisre;
 260     }
 261     ndo_real(drblock->rav,drblock->ndr,bOK1);
 262     bOK = bOK && bOK1;
 263     ndo_real(drblock->rt,drblock->ndr,bOK1);
 264     bOK = bOK && bOK1;
 265     if (bRead)
 266       *ndr = eh.ndisre;
 267   }
 268   else
 269     if (bRead)
 270       *ndr = 0;
 271
 272   if (eh.u_size) {
 273     fatal_error(0,"Can't handle user blocks");
 274   }
 275   if (!bOK) {
 276     if (bRead) {
 277       fprintf(stderr,"\nLast frame read %d                               ",framenr-1);
 278       fprintf(stderr,"\nWARNING: Incomplete frame: nr %6d time %8.3f     \n",
 279               framenr,*t);
 280     } else
 281       fatal_error(-1,"could not write energies");
 282     return FALSE;
 283   }
 284
 285   return TRUE;
 286 }
 287
 288
 289