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 3.1
  11  * Copyright (c) 1991-2001, University of Groningen, The Netherlands
  12  * This program is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU General Public License
  14  * as published by the Free Software Foundation; either version 2
  15  * of the License, or (at your option) any later version.
  16  *
  17  * If you want to redistribute modifications, please consider that
  18  * scientific software is very special. Version control is crucial -
  19  * bugs must be traceable. We will be happy to consider code for
  20  * inclusion in the official distribution, but derived work must not
  21  * be called official GROMACS. Details are found in the README & COPYING
  22  * files - if they are missing, get the official version at www.gromacs.org.
  23  *
  24  * To help us fund GROMACS development, we humbly ask that you cite
  25  * the papers on the package - you can find them in the top README file.
  26  *
  27  * For more info, check our website at http://www.gromacs.org
  28  *
  29  * And Hey:
  30  * Great Red Owns Many ACres of Sand
  31  */
  32
  33 #include "futil.h"
  34 #include "string2.h"
  35 #include "fatal.h"
  36 #include "smalloc.h"
  37 #include "gmxfio.h"
  38 #include "enxio.h"
  39
  40 void free_enxframe(t_enxframe *fr)
  41 {
  42   int b;
  43
  44   if (fr->e_alloc)
  45     sfree(fr->ener);
  46   if (fr->d_alloc) {
  47     sfree(fr->rav);
  48     sfree(fr->rt);
  49   }
  50   for(b=0; b<fr->nblock; b++)
  51     sfree(fr->block[b]);
  52   sfree(fr->block);
  53   sfree(fr->b_alloc);
  54   sfree(fr->nr);
  55 }
  56
  57 static void wr_ener_nms(FILE *out,int nre,char *nms[])
  58 {
  59   int i;
  60
  61   fprintf(out,"%d\n",nre);
  62   for(i=0; (i<nre); i++)
  63     fprintf(out,"%s\n",nms[i]);
  64 }
  65
  66 static void rd_ener_nms(FILE *in,int *nre,char ***nm)
  67 {
  68   char line[256];
  69   int  i;
  70
  71   fgets2(line,255,in);
  72   if (sscanf(line,"%d",nre) == 0) {
  73     *nre=0;
  74     return;
  75   }
  76   snew((*nm),*nre);
  77   for(i=0; (i< (*nre)); i++) {
  78     fgets2(line,255,in);
  79     trim(line);
  80     (*nm)[i]=strdup(line);
  81   }
  82 }
  83
  84 static void edr_nms(int fp,int *nre,char ***nms)
  85 {
  86   XDR  *xdr;
  87   bool bRead = fio_getread(fp);
  88   int  i;
  89   char **NM;
  90
  91   xdr = fio_getxdr(fp);
  92
  93   NM=*nms;
  94
  95   if (!xdr_int(xdr,nre)) {
  96     *nre=0;
  97     return;
  98   }
  99   if (NM == NULL) {
 100     snew(NM,*nre);
 101   }
 102   for(i=0; i<*nre; i++) {
 103     if (bRead && NM[i]) {
 104       sfree(NM[i]);
 105       NM[i] = NULL;
 106     }
 107     xdr_string(xdr,&(NM[i]),STRLEN);
 108   }
 109   *nms=NM;
 110 }
 111
 112 static bool do_eheader(int fp,t_enxframe *fr,bool *bOK)
 113 {
 114   int  block,i,dum=0;
 115   bool bRead = fio_getread(fp);
 116   int  tempfix_nr=0;
 117
 118   *bOK=TRUE;
 119   if (!do_real(fr->t))       return FALSE;
 120   if (!do_int (fr->step))    *bOK = FALSE;
 121   if (!do_int (fr->nre))     *bOK = FALSE;
 122   if (!do_int (fr->ndisre))  *bOK = FALSE;
 123   if (!do_int (fr->nblock))  *bOK = FALSE;
 124   if (bRead && fr->nblock>70) {
 125     /* Temporary fix for intermediate file version, only used by B. Hess */
 126     tempfix_nr = fr->nblock;
 127     fr->nblock = 1;
 128   }
 129   if (bRead && fr->nblock>fr->nr_alloc) {
 130     srenew(fr->nr,fr->nblock);
 131     srenew(fr->b_alloc,fr->nblock);
 132     srenew(fr->block,fr->nblock);
 133     for(i=fr->nr_alloc; i<fr->nblock; i++) {
 134       fr->block[i]   = NULL;
 135       fr->b_alloc[i] = 0;
 136     }
 137     fr->nr_alloc = fr->nblock;
 138   }
 139   if (tempfix_nr)
 140     fr->nr[0] = tempfix_nr;
 141   else {
 142     for(block=0; block<fr->nblock; block++)
 143       if (!do_int (fr->nr[block]))
 144         *bOK = FALSE;
 145   }
 146   if (!do_int (fr->e_size))  *bOK = FALSE;
 147   if (!do_int (fr->d_size))  *bOK = FALSE;
 148   /* Do a dummy int to keep the format compatible with the old code */
 149   if (!do_int (dum))         *bOK = FALSE;
 150
 151   return *bOK;
 152 }
 153
 154 void do_enxnms(int fp,int *nre,char ***nms)
 155 {
 156   bool bRead;
 157
 158   bRead = fio_getread(fp);
 159   if (fio_getftp(fp) == efEDR) {
 160     fio_select(fp);
 161     edr_nms(fp,nre,nms);
 162   }
 163   else if (bRead)
 164     rd_ener_nms(fio_getfp(fp),nre,nms);
 165   else
 166     wr_ener_nms(fio_getfp(fp),*nre,*nms);
 167 }
 168
 169 void close_enx(int fp)
 170 {
 171   fio_close(fp);
 172 }
 173
 174 static bool empty_file(char *fn)
 175 {
 176   FILE *fp;
 177   char dum;
 178   bool bEmpty;
 179
 180   fp = ffopen(fn,"r");
 181   fread(&dum,sizeof(dum),1,fp);
 182   bEmpty = feof(fp);
 183   fclose(fp);
 184
 185   return bEmpty;
 186 }
 187
 188 static int framenr;
 189
 190 int open_enx(char *fn,char *mode)
 191 {
 192   int        fp,nre,i;
 193   char       **nm=NULL;
 194   t_enxframe *fr;
 195   bool       bDum=TRUE;
 196
 197   if (mode[0]=='r') {
 198     fp=fio_open(fn,mode);
 199     fio_select(fp);
 200     fio_setprecision(fp,FALSE);
 201     do_enxnms(fp,&nre,&nm);
 202     snew(fr,1);
 203     do_eheader(fp,fr,&bDum);
 204
 205     /* Now check whether this file is in single precision */
 206     if (((fr->e_size && (fr->nre == nre) &&
 207           (nre*4*sizeof(float) == fr->e_size)) ||
 208          (fr->d_size &&
 209           (fr->ndisre*sizeof(float)*2+sizeof(int) == fr->d_size)))){
 210       fprintf(stderr,"Opened %s as single precision energy file\n",fn);
 211       for(i=0; (i<nre); i++)
 212         sfree(nm[i]);
 213       sfree(nm);
 214     }
 215     else {
 216       fio_rewind(fp);
 217       fio_select(fp);
 218       fio_setprecision(fp,TRUE);
 219       do_enxnms(fp,&nre,&nm);
 220       do_eheader(fp,fr,&bDum);
 221       if (((fr->e_size && (fr->nre == nre) &&
 222             (nre*4*sizeof(double) == fr->e_size)) ||
 223            (fr->d_size &&
 224             (fr->ndisre*sizeof(double)*2+sizeof(int) == fr->d_size))))
 225         fprintf(stderr,"Opened %s as double precision energy file\n",fn);
 226       else {
 227         if (empty_file(fn))
 228           fatal_error(0,"File %s is empty",fn);
 229         else
 230           fatal_error(0,"Energy file %s not recognized, maybe different CPU?",
 231                       fn);
 232       }
 233       for(i=0; (i<nre); i++)
 234           sfree(nm[i]);
 235       sfree(nm);
 236     }
 237     free_enxframe(fr);
 238     sfree(fr);
 239     fio_rewind(fp);
 240   }
 241   else
 242     fp = fio_open(fn,mode);
 243
 244   framenr=0;
 245
 246   return fp;
 247 }
 248
 249 bool do_enx(int fp,t_enxframe *fr)
 250 {
 251   int       i,block;
 252   bool      bRead,bOK,bOK1,bSane;
 253   real      tmp1,tmp2;
 254
 255   bOK = TRUE;
 256   bRead = fio_getread(fp);
 257   if (!bRead) {
 258     fr->e_size = fr->nre*sizeof(fr->ener[0].e)*4;
 259     fr->d_size = fr->ndisre*(sizeof(fr->rav[0]) +
 260                            sizeof(fr->rt[0]));
 261   }
 262   fio_select(fp);
 263
 264   if (!do_eheader(fp,fr,&bOK)) {
 265     if (bRead) {
 266         fprintf(stderr,"\rLast frame read %d                          ",framenr-1);
 267         if (!bOK)
 268           fprintf(stderr,"\nWARNING: Incomplete frame: nr %6d time %8.3f\n",
 269                   framenr,fr->t);
 270     }
 271     return FALSE;
 272   }
 273   if (bRead) {
 274     if ( ( framenr<10 ) || ( framenr%10 == 0) )
 275       fprintf(stderr,"\rReading frame %6d time %8.3f           ",framenr,fr->t);
 276     framenr++;
 277   }
 278   /* Check sanity of this header */
 279   bSane = (fr->nre > 0 || fr->ndisre > 0);
 280   for(block=0; block<fr->nblock; block++)
 281     bSane = bSane || (fr->nr[block] > 0);
 282   if (!((fr->step >= 0) && bSane)) {
 283     fprintf(stderr,"\nWARNING: there may be something wrong with energy file %s\n",
 284             fio_getname(fp));
 285     fprintf(stderr,"Found: step=%d, nre=%d, ndisre=%d, nblock=%d, time=%g.\n"
 286             "Trying to skip frame expect a crash though\n",
 287             fr->step,fr->nre,fr->ndisre,fr->nblock,fr->t);
 288   }
 289   if (bRead && fr->nre>fr->e_alloc) {
 290     srenew(fr->ener,fr->nre);
 291     fr->e_alloc = fr->nre;
 292   }
 293   for(i=0; i<fr->nre; i++) {
 294     bOK = bOK && do_real(fr->ener[i].e);
 295
 296     tmp1 = fr->ener[i].eav;
 297     if((tmp1/(fr->step+1))<GMX_REAL_EPS)
 298       tmp1=0;
 299     bOK = bOK && do_real(tmp1);
 300     fr->ener[i].eav = tmp1;
 301
 302     /* This is to save only in single precision (unless compiled in DP) */
 303     tmp2 = fr->ener[i].esum;
 304     bOK = bOK && do_real(tmp2);
 305     fr->ener[i].esum = tmp2;
 306
 307     bOK = bOK && do_real(fr->ener[i].e2sum);
 308   }
 309   if (fr->ndisre) {
 310     if (bRead && fr->ndisre>fr->d_alloc) {
 311       srenew(fr->rav,fr->ndisre);
 312       srenew(fr->rt,fr->ndisre);
 313       fr->d_alloc = fr->ndisre;
 314     }
 315     ndo_real(fr->rav,fr->ndisre,bOK1);
 316     bOK = bOK && bOK1;
 317     ndo_real(fr->rt,fr->ndisre,bOK1);
 318     bOK = bOK && bOK1;
 319   }
 320   for(block=0; block<fr->nblock; block++) {
 321     if (bRead && fr->nr[block]>fr->b_alloc[block]) {
 322       srenew(fr->block[block],fr->nr[block]);
 323       fr->b_alloc[block] = fr->nr[block];
 324     }
 325     ndo_real(fr->block[block],fr->nr[block],bOK1);
 326     bOK = bOK && bOK1;
 327   }
 328   if (!bOK) {
 329     if (bRead) {
 330       fprintf(stderr,"\nLast frame read %d                               ",
 331               framenr-1);
 332       fprintf(stderr,"\nWARNING: Incomplete frame: nr %6d time %8.3f     \n",
 333               framenr,fr->t);
 334     } else
 335       fatal_error(-1,"could not write energies");
 336     return FALSE;
 337   }
 338
 339   return TRUE;
 340 }
 341
 342
 343