src/gmxlib/enxio.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  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include "futil.h"
  43 #include "string2.h"
  44 #include "gmx_fatal.h"
  45 #include "smalloc.h"
  46 #include "gmxfio.h"
  47 #include "enxio.h"
  48 #include "vec.h"
  49 #include "xdrf.h"
  50
  51 /* The source code in this file should be thread-safe.
  52          Please keep it that way. */
  53
  54 /* This number should be increased whenever the file format changes! */
  55 static const int enx_version = 5;
  56
  57 const char *enx_block_id_name[] = {
  58     "Averaged orientation restraints",
  59     "Instantaneous orientation restraints",
  60     "Orientation restraint order tensor(s)",
  61     "Distance restraints",
  62     "Free energy data",
  63     "BAR histogram",
  64     "Delta H raw data"
  65 };
  66
  67
  68 /* Stuff for reading pre 4.1 energy files */
  69 typedef struct {
  70     gmx_bool     bOldFileOpen;   /* Is this an open old file? */
  71     gmx_bool     bReadFirstStep; /* Did we read the first step? */
  72     int      first_step;     /* First step in the energy file */
  73     int      step_prev;      /* Previous step */
  74     int      nsum_prev;      /* Previous step sum length */
  75     t_energy *ener_prev;     /* Previous energy sums */
  76 } ener_old_t;
  77
  78 struct ener_file
  79 {
  80     ener_old_t eo;
  81     t_fileio *fio;
  82     int framenr;
  83     real frametime;
  84 };
  85
  86 static void enxsubblock_init(t_enxsubblock *sb)
  87 {
  88     sb->nr=0;
  89 #ifdef GMX_DOUBLE
  90     sb->type=xdr_datatype_double;
  91 #else
  92     sb->type=xdr_datatype_float;
  93 #endif
  94     sb->fval = NULL;
  95     sb->dval = NULL;
  96     sb->ival = NULL;
  97     sb->lval = NULL;
  98     sb->cval = NULL;
  99     sb->sval = NULL;
 100     sb->fval_alloc = 0;
 101     sb->dval_alloc = 0;
 102     sb->ival_alloc = 0;
 103     sb->lval_alloc = 0;
 104     sb->cval_alloc = 0;
 105     sb->sval_alloc = 0;
 106 }
 107
 108 static void enxsubblock_free(t_enxsubblock *sb)
 109 {
 110     if (sb->fval_alloc)
 111     {
 112         free(sb->fval);
 113         sb->fval_alloc=0;
 114         sb->fval=NULL;
 115     }
 116     if (sb->dval_alloc)
 117     {
 118         free(sb->dval);
 119         sb->dval_alloc=0;
 120         sb->dval=NULL;
 121     }
 122     if (sb->ival_alloc)
 123     {
 124         free(sb->ival);
 125         sb->ival_alloc=0;
 126         sb->ival=NULL;
 127     }
 128     if (sb->lval_alloc)
 129     {
 130         free(sb->lval);
 131         sb->lval_alloc=0;
 132         sb->lval=NULL;
 133     }
 134     if (sb->cval_alloc)
 135     {
 136         free(sb->cval);
 137         sb->cval_alloc=0;
 138         sb->cval=NULL;
 139     }
 140     if (sb->sval_alloc)
 141     {
 142         int i;
 143
 144         for(i=0;i<sb->sval_alloc;i++)
 145         {
 146             if (sb->sval[i])
 147             {
 148                 free(sb->sval[i]);
 149             }
 150         }
 151         free(sb->sval);
 152         sb->sval_alloc=0;
 153         sb->sval=NULL;
 154     }
 155 }
 156
 157 /* allocate the appropriate amount of memory for the given type and nr */
 158 static void enxsubblock_alloc(t_enxsubblock *sb)
 159 {
 160     /* allocate the appropriate amount of memory */
 161     switch(sb->type)
 162     {
 163         case xdr_datatype_float:
 164             if (sb->nr > sb->fval_alloc)
 165             {
 166                 srenew(sb->fval, sb->nr);
 167                 sb->fval_alloc=sb->nr;
 168             }
 169             break;
 170         case xdr_datatype_double:
 171             if (sb->nr > sb->dval_alloc)
 172             {
 173                 srenew(sb->dval, sb->nr);
 174                 sb->dval_alloc=sb->nr;
 175             }
 176             break;
 177         case xdr_datatype_int:
 178             if (sb->nr > sb->ival_alloc)
 179             {
 180                 srenew(sb->ival, sb->nr);
 181                 sb->ival_alloc=sb->nr;
 182             }
 183             break;
 184         case xdr_datatype_large_int:
 185             if (sb->nr > sb->lval_alloc)
 186             {
 187                 srenew(sb->lval, sb->nr);
 188                 sb->lval_alloc=sb->nr;
 189             }
 190             break;
 191         case xdr_datatype_char:
 192             if (sb->nr > sb->cval_alloc)
 193             {
 194                 srenew(sb->cval, sb->nr);
 195                 sb->cval_alloc=sb->nr;
 196             }
 197             break;
 198         case xdr_datatype_string:
 199             if (sb->nr > sb->sval_alloc)
 200             {
 201                 int i;
 202
 203                 srenew(sb->sval, sb->nr);
 204                 for(i=sb->sval_alloc;i<sb->nr;i++)
 205                 {
 206                     sb->sval[i]=NULL;
 207                 }
 208                 sb->sval_alloc=sb->nr;
 209             }
 210             break;
 211         default:
 212             gmx_incons("Unknown block type: this file is corrupted or from the future");
 213     }
 214 }
 215
 216 static void enxblock_init(t_enxblock *eb)
 217 {
 218     eb->id=enxOR;
 219     eb->nsub=0;
 220     eb->sub=NULL;
 221     eb->nsub_alloc=0;
 222 }
 223
 224 static void enxblock_free(t_enxblock *eb)
 225 {
 226     if (eb->nsub_alloc>0)
 227     {
 228         int i;
 229         for(i=0;i<eb->nsub_alloc;i++)
 230         {
 231             enxsubblock_free(&(eb->sub[i]));
 232         }
 233         free(eb->sub);
 234         eb->nsub_alloc=0;
 235         eb->sub=NULL;
 236     }
 237 }
 238
 239 void init_enxframe(t_enxframe *fr)
 240 {
 241     fr->e_alloc=0;
 242     fr->ener=NULL;
 243
 244     /*fr->d_alloc=0;*/
 245     fr->ener=NULL;
 246
 247     /*fr->ndisre=0;*/
 248
 249     fr->nblock=0;
 250     fr->nblock_alloc=0;
 251     fr->block=NULL;
 252 }
 253
 254
 255 void free_enxframe(t_enxframe *fr)
 256 {
 257   int b;
 258
 259   if (fr->e_alloc)
 260   {
 261     sfree(fr->ener);
 262   }
 263   for(b=0; b<fr->nblock_alloc; b++)
 264   {
 265       enxblock_free(&(fr->block[b]));
 266   }
 267   free(fr->block);
 268 }
 269
 270 void add_blocks_enxframe(t_enxframe *fr, int n)
 271 {
 272     fr->nblock=n;
 273     if (n > fr->nblock_alloc)
 274     {
 275         int b;
 276
 277         srenew(fr->block, n);
 278         for(b=fr->nblock_alloc;b<fr->nblock;b++)
 279         {
 280             enxblock_init(&(fr->block[b]));
 281         }
 282         fr->nblock_alloc=n;
 283     }
 284 }
 285
 286 t_enxblock *find_block_id_enxframe(t_enxframe *ef, int id, t_enxblock *prev)
 287 {
 288     gmx_off_t starti=0;
 289     gmx_off_t i;
 290
 291     if (prev)
 292     {
 293         starti=(prev - ef->block) + 1;
 294     }
 295     for(i=starti; i<ef->nblock; i++)
 296     {
 297         if (ef->block[i].id == id)
 298             return &(ef->block[i]);
 299     }
 300     return NULL;
 301 }
 302
 303 void add_subblocks_enxblock(t_enxblock *eb, int n)
 304 {
 305     eb->nsub=n;
 306     if (eb->nsub > eb->nsub_alloc)
 307     {
 308         int b;
 309
 310         srenew(eb->sub, n);
 311         for(b=eb->nsub_alloc; b<n; b++)
 312         {
 313             enxsubblock_init(&(eb->sub[b]));
 314         }
 315         eb->nsub_alloc=n;
 316     }
 317 }
 318
 319 static void enx_warning(const char *msg)
 320 {
 321     if (getenv("GMX_ENX_NO_FATAL") != NULL)
 322     {
 323         gmx_warning(msg);
 324     }
 325     else
 326     {
 327         gmx_fatal(FARGS,"%s\n%s",
 328                   msg,
 329                   "If you want to use the correct frames before the corrupted frame and avoid this fatal error set the env.var. GMX_ENX_NO_FATAL");
 330     }
 331 }
 332
 333 static void edr_strings(XDR *xdr,gmx_bool bRead,int file_version,
 334                         int n,gmx_enxnm_t **nms)
 335 {
 336     int  i;
 337     gmx_enxnm_t *nm;
 338
 339     if (*nms == NULL)
 340     {
 341         snew(*nms,n);
 342     }
 343     for(i=0; i<n; i++)
 344     {
 345         nm = &(*nms)[i];
 346         if (bRead)
 347         {
 348             if (nm->name)
 349             {
 350                 sfree(nm->name);
 351                 nm->name = NULL;
 352             }
 353             if (nm->unit)
 354             {
 355                 sfree(nm->unit);
 356                 nm->unit = NULL;
 357             }
 358         }
 359         if(!xdr_string(xdr,&(nm->name),STRLEN))
 360         {
 361             gmx_file("Cannot write energy names to file; maybe you are out of disk space?");
 362         }
 363         if (file_version >= 2)
 364         {
 365             if(!xdr_string(xdr,&(nm->unit),STRLEN))
 366             {
 367                 gmx_file("Cannot write energy names to file; maybe you are out of disk space?");
 368             }
 369         }
 370         else
 371         {
 372             nm->unit = strdup("kJ/mol");
 373         }
 374     }
 375 }
 376
 377 void do_enxnms(ener_file_t ef,int *nre,gmx_enxnm_t **nms)
 378 {
 379     int  magic=-55555;
 380     XDR  *xdr;
 381     gmx_bool bRead = gmx_fio_getread(ef->fio);
 382     int  file_version;
 383     int  i;
 384
 385     gmx_fio_checktype(ef->fio);
 386
 387     xdr = gmx_fio_getxdr(ef->fio);
 388
 389     if (!xdr_int(xdr,&magic))
 390     {
 391         if(!bRead)
 392         {
 393             gmx_file("Cannot write energy names to file; maybe you are out of disk space?");
 394         }
 395         *nre=0;
 396         return;
 397     }
 398     if (magic > 0)
 399     {
 400         /* Assume this is an old edr format */
 401         file_version = 1;
 402         *nre = magic;
 403         ef->eo.bOldFileOpen = TRUE;
 404         ef->eo.bReadFirstStep = FALSE;
 405         srenew(ef->eo.ener_prev,*nre);
 406     }
 407     else
 408     {
 409         ef->eo.bOldFileOpen=FALSE;
 410
 411         if (magic != -55555)
 412         {
 413             gmx_fatal(FARGS,"Energy names magic number mismatch, this is not a GROMACS edr file");
 414         }
 415         file_version = enx_version;
 416         xdr_int(xdr,&file_version);
 417         if (file_version > enx_version)
 418         {
 419             gmx_fatal(FARGS,"reading tpx file (%s) version %d with version %d program",gmx_fio_getname(ef->fio),file_version,enx_version);
 420         }
 421         xdr_int(xdr,nre);
 422     }
 423     if (file_version != enx_version)
 424     {
 425         fprintf(stderr,"Note: enx file_version %d, software version %d\n",
 426                 file_version,enx_version);
 427     }
 428
 429     edr_strings(xdr,bRead,file_version,*nre,nms);
 430 }
 431
 432 static gmx_bool do_eheader(ener_file_t ef,int *file_version,t_enxframe *fr,
 433                        int nre_test,gmx_bool *bWrongPrecision,gmx_bool *bOK)
 434 {
 435     int  magic=-7777777;
 436     real first_real_to_check;
 437     int  b,i,zero=0,dum=0;
 438     gmx_bool bRead = gmx_fio_getread(ef->fio);
 439     int  tempfix_nr=0;
 440     int  ndisre=0;
 441     int  startb=0;
 442 #ifndef GMX_DOUBLE
 443     xdr_datatype dtreal=xdr_datatype_float;
 444 #else
 445     xdr_datatype dtreal=xdr_datatype_double;
 446 #endif
 447
 448     if (bWrongPrecision)
 449     {
 450         *bWrongPrecision = FALSE;
 451     }
 452
 453     *bOK=TRUE;
 454     /* The original energy frame started with a real,
 455      * so we have to use a real for compatibility.
 456      * This is VERY DIRTY code, since do_eheader can be called
 457      * with the wrong precision set and then we could read r > -1e10,
 458      * while actually the intention was r < -1e10.
 459      * When nre_test >= 0, do_eheader should therefore terminate
 460      * before the number of i/o calls starts depending on what has been read
 461      * (which is the case for for instance the block sizes for variable
 462      * number of blocks, where this number is read before).
 463      */
 464     first_real_to_check = -2e10;
 465     if (!gmx_fio_do_real(ef->fio, first_real_to_check))
 466     {
 467         return FALSE;
 468     }
 469     if (first_real_to_check > -1e10)
 470     {
 471         /* Assume we are reading an old format */
 472         *file_version = 1;
 473         fr->t = first_real_to_check;
 474         if (!gmx_fio_do_int(ef->fio, dum))   *bOK = FALSE;
 475         fr->step = dum;
 476     }
 477     else
 478     {
 479         if (!gmx_fio_do_int(ef->fio, magic))       *bOK = FALSE;
 480         if (magic != -7777777)
 481         {
 482             enx_warning("Energy header magic number mismatch, this is not a GROMACS edr file");
 483             *bOK=FALSE;
 484             return FALSE;
 485         }
 486         *file_version = enx_version;
 487         if (!gmx_fio_do_int(ef->fio, *file_version)) *bOK = FALSE;
 488         if (*bOK && *file_version > enx_version)
 489         {
 490             gmx_fatal(FARGS,"reading tpx file (%s) version %d with version %d program",gmx_fio_getname(ef->fio),file_version,enx_version);
 491         }
 492         if (!gmx_fio_do_double(ef->fio, fr->t))       *bOK = FALSE;
 493         if (!gmx_fio_do_gmx_large_int(ef->fio, fr->step)) *bOK = FALSE;
 494         if (!bRead && fr->nsum == 1) {
 495             /* Do not store sums of length 1,
 496              * since this does not add information.
 497              */
 498             if (!gmx_fio_do_int(ef->fio, zero))      *bOK = FALSE;
 499         } else {
 500             if (!gmx_fio_do_int(ef->fio, fr->nsum))  *bOK = FALSE;
 501         }
 502         if (*file_version >= 3)
 503         {
 504             if (!gmx_fio_do_gmx_large_int(ef->fio, fr->nsteps)) *bOK = FALSE;
 505         }
 506         else
 507         {
 508             fr->nsteps = max(1,fr->nsum);
 509         }
 510         if (*file_version >= 5)
 511         {
 512             if (!gmx_fio_do_double(ef->fio, fr->dt)) *bOK = FALSE;
 513         }
 514         else
 515         {
 516             fr->dt = 0;
 517         }
 518     }
 519     if (!gmx_fio_do_int(ef->fio, fr->nre))     *bOK = FALSE;
 520     if (*file_version < 4)
 521     {
 522         if (!gmx_fio_do_int(ef->fio, ndisre))  *bOK = FALSE;
 523     }
 524     else
 525     {
 526         /* now reserved for possible future use */
 527         if (!gmx_fio_do_int(ef->fio, dum))  *bOK = FALSE;
 528     }
 529
 530     if (!gmx_fio_do_int(ef->fio, fr->nblock))  *bOK = FALSE;
 531     if (fr->nblock < 0) *bOK=FALSE;
 532
 533     if (ndisre!=0)
 534     {
 535         if (*file_version >= 4)
 536         {
 537             enx_warning("Distance restraint blocks in old style in new style file");
 538             *bOK=FALSE;
 539             return FALSE;
 540         }
 541         fr->nblock+=1;
 542     }
 543
 544
 545     /* Frames could have nre=0, so we can not rely only on the fr->nre check */
 546     if (bRead && nre_test >= 0 &&
 547         ((fr->nre > 0 && fr->nre != nre_test) ||
 548          fr->nre < 0 || ndisre < 0 || fr->nblock < 0))
 549     {
 550         *bWrongPrecision = TRUE;
 551         return *bOK;
 552     }
 553
 554     /* we now know what these should be, or we've already bailed out because
 555        of wrong precision */
 556     if ( *file_version==1 && (fr->t < 0 || fr->t > 1e20 || fr->step < 0 ) )
 557     {
 558         enx_warning("edr file with negative step number or unreasonable time (and without version number).");
 559         *bOK=FALSE;
 560         return FALSE;
 561     }
 562
 563
 564     if (*bOK && bRead)
 565     {
 566         add_blocks_enxframe(fr, fr->nblock);
 567     }
 568
 569     startb=0;
 570     if (ndisre>0)
 571     {
 572         /* sub[0] is the instantaneous data, sub[1] is time averaged */
 573         add_subblocks_enxblock(&(fr->block[0]), 2);
 574         fr->block[0].id=enxDISRE;
 575         fr->block[0].sub[0].nr=ndisre;
 576         fr->block[0].sub[1].nr=ndisre;
 577         fr->block[0].sub[0].type=dtreal;
 578         fr->block[0].sub[1].type=dtreal;
 579         startb++;
 580     }
 581
 582     /* read block header info */
 583     for(b=startb; b<fr->nblock; b++)
 584     {
 585         if (*file_version<4)
 586         {
 587             /* blocks in old version files always have 1 subblock that
 588                consists of reals. */
 589             int nrint;
 590
 591             if (bRead)
 592             {
 593                 add_subblocks_enxblock(&(fr->block[b]), 1);
 594             }
 595             else
 596             {
 597                 if (fr->block[b].nsub != 1)
 598                 {
 599                     gmx_incons("Writing an old version .edr file with too many subblocks");
 600                 }
 601                 if (fr->block[b].sub[0].type != dtreal)
 602                 {
 603                     gmx_incons("Writing an old version .edr file the wrong subblock type");
 604                 }
 605             }
 606             nrint = fr->block[b].sub[0].nr;
 607
 608             if (!gmx_fio_do_int(ef->fio, nrint))
 609             {
 610                 *bOK = FALSE;
 611             }
 612             fr->block[b].id          = b - startb;
 613             fr->block[b].sub[0].nr   = nrint;
 614             fr->block[b].sub[0].type = dtreal;
 615         }
 616         else
 617         {
 618             int i;
 619             /* in the new version files, the block header only contains
 620                the ID and the number of subblocks */
 621             int nsub=fr->block[b].nsub;
 622             *bOK = *bOK && gmx_fio_do_int(ef->fio, fr->block[b].id);
 623             *bOK = *bOK && gmx_fio_do_int(ef->fio, nsub);
 624
 625             fr->block[b].nsub=nsub;
 626             if (bRead)
 627                 add_subblocks_enxblock(&(fr->block[b]), nsub);
 628
 629             /* read/write type & size for each subblock */
 630             for(i=0;i<nsub;i++)
 631             {
 632                 t_enxsubblock *sub=&(fr->block[b].sub[i]); /* shortcut */
 633                 int typenr=sub->type;
 634
 635                 *bOK=*bOK && gmx_fio_do_int(ef->fio, typenr);
 636                 *bOK=*bOK && gmx_fio_do_int(ef->fio, sub->nr);
 637
 638                 sub->type = (xdr_datatype)typenr;
 639             }
 640         }
 641     }
 642     if (!gmx_fio_do_int(ef->fio, fr->e_size))  *bOK = FALSE;
 643
 644     /* now reserved for possible future use */
 645     if (!gmx_fio_do_int(ef->fio, dum))  *bOK = FALSE;
 646
 647     /* Do a dummy int to keep the format compatible with the old code */
 648     if (!gmx_fio_do_int(ef->fio, dum))         *bOK = FALSE;
 649
 650     if (*bOK && *file_version == 1 && nre_test < 0)
 651     {
 652 #if 0
 653         if (fp >= ener_old_nalloc)
 654         {
 655             gmx_incons("Problem with reading old format energy files");
 656         }
 657 #endif
 658
 659         if (!ef->eo.bReadFirstStep)
 660         {
 661             ef->eo.bReadFirstStep = TRUE;
 662             ef->eo.first_step     = fr->step;
 663             ef->eo.step_prev      = fr->step;
 664             ef->eo.nsum_prev      = 0;
 665         }
 666
 667         fr->nsum   = fr->step - ef->eo.first_step + 1;
 668         fr->nsteps = fr->step - ef->eo.step_prev;
 669         fr->dt     = 0;
 670     }
 671
 672     return *bOK;
 673 }
 674
 675 void free_enxnms(int n,gmx_enxnm_t *nms)
 676 {
 677     int i;
 678
 679     for(i=0; i<n; i++)
 680     {
 681         sfree(nms[i].name);
 682         sfree(nms[i].unit);
 683     }
 684
 685     sfree(nms);
 686 }
 687
 688 void close_enx(ener_file_t ef)
 689 {
 690     if(gmx_fio_close(ef->fio) != 0)
 691     {
 692         gmx_file("Cannot close energy file; it might be corrupt, or maybe you are out of disk space?");
 693     }
 694 }
 695
 696 static gmx_bool empty_file(const char *fn)
 697 {
 698     FILE *fp;
 699     char dum;
 700     int  ret;
 701     gmx_bool bEmpty;
 702
 703     fp = gmx_fio_fopen(fn,"r");
 704     ret = fread(&dum,sizeof(dum),1,fp);
 705     bEmpty = feof(fp);
 706     gmx_fio_fclose(fp);
 707
 708     return bEmpty;
 709 }
 710
 711
 712 ener_file_t open_enx(const char *fn,const char *mode)
 713 {
 714     int        nre,i;
 715     gmx_enxnm_t *nms=NULL;
 716     int        file_version=-1;
 717     t_enxframe *fr;
 718     gmx_bool       bWrongPrecision,bOK=TRUE;
 719     struct ener_file *ef;
 720
 721     snew(ef,1);
 722
 723     if (mode[0]=='r') {
 724         ef->fio=gmx_fio_open(fn,mode);
 725         gmx_fio_checktype(ef->fio);
 726         gmx_fio_setprecision(ef->fio,FALSE);
 727         do_enxnms(ef,&nre,&nms);
 728         snew(fr,1);
 729         do_eheader(ef,&file_version,fr,nre,&bWrongPrecision,&bOK);
 730         if(!bOK)
 731         {
 732             gmx_file("Cannot read energy file header. Corrupt file?");
 733         }
 734
 735         /* Now check whether this file is in single precision */
 736         if (!bWrongPrecision &&
 737             ((fr->e_size && (fr->nre == nre) &&
 738               (nre*4*(long int)sizeof(float) == fr->e_size)) ) )
 739         {
 740             fprintf(stderr,"Opened %s as single precision energy file\n",fn);
 741             free_enxnms(nre,nms);
 742         }
 743         else
 744         {
 745             gmx_fio_rewind(ef->fio);
 746             gmx_fio_checktype(ef->fio);
 747             gmx_fio_setprecision(ef->fio,TRUE);
 748             do_enxnms(ef,&nre,&nms);
 749             do_eheader(ef,&file_version,fr,nre,&bWrongPrecision,&bOK);
 750             if(!bOK)
 751             {
 752                 gmx_file("Cannot write energy file header; maybe you are out of disk space?");
 753             }
 754
 755             if (((fr->e_size && (fr->nre == nre) &&
 756                             (nre*4*(long int)sizeof(double) == fr->e_size)) ))
 757                 fprintf(stderr,"Opened %s as double precision energy file\n",
 758                         fn);
 759             else {
 760                 if (empty_file(fn))
 761                     gmx_fatal(FARGS,"File %s is empty",fn);
 762                 else
 763                     gmx_fatal(FARGS,"Energy file %s not recognized, maybe different CPU?",
 764                               fn);
 765             }
 766             free_enxnms(nre,nms);
 767         }
 768         free_enxframe(fr);
 769         sfree(fr);
 770         gmx_fio_rewind(ef->fio);
 771     }
 772     else
 773         ef->fio = gmx_fio_open(fn,mode);
 774
 775     ef->framenr=0;
 776     ef->frametime=0;
 777     return ef;
 778 }
 779
 780 t_fileio *enx_file_pointer(const ener_file_t ef)
 781 {
 782     return ef->fio;
 783 }
 784
 785 static void convert_full_sums(ener_old_t *ener_old,t_enxframe *fr)
 786 {
 787     int nstep_all;
 788     int ne,ns,i;
 789     double esum_all,eav_all;
 790
 791     if (fr->nsum > 0)
 792     {
 793         ne = 0;
 794         ns = 0;
 795         for(i=0; i<fr->nre; i++)
 796         {
 797             if (fr->ener[i].e    != 0) ne++;
 798             if (fr->ener[i].esum != 0) ns++;
 799         }
 800         if (ne > 0 && ns == 0)
 801         {
 802             /* We do not have all energy sums */
 803             fr->nsum = 0;
 804         }
 805     }
 806
 807     /* Convert old full simulation sums to sums between energy frames */
 808     nstep_all = fr->step - ener_old->first_step + 1;
 809     if (fr->nsum > 1 && fr->nsum == nstep_all && ener_old->nsum_prev > 0)
 810     {
 811         /* Set the new sum length: the frame step difference */
 812         fr->nsum = fr->step - ener_old->step_prev;
 813         for(i=0; i<fr->nre; i++)
 814         {
 815             esum_all = fr->ener[i].esum;
 816             eav_all  = fr->ener[i].eav;
 817             fr->ener[i].esum = esum_all - ener_old->ener_prev[i].esum;
 818             fr->ener[i].eav  = eav_all  - ener_old->ener_prev[i].eav
 819                 - dsqr(ener_old->ener_prev[i].esum/(nstep_all - fr->nsum)
 820                        - esum_all/nstep_all)*
 821                 (nstep_all - fr->nsum)*nstep_all/(double)fr->nsum;
 822             ener_old->ener_prev[i].esum = esum_all;
 823             ener_old->ener_prev[i].eav  = eav_all;
 824         }
 825         ener_old->nsum_prev = nstep_all;
 826     }
 827     else if (fr->nsum > 0)
 828     {
 829         if (fr->nsum != nstep_all)
 830         {
 831             fprintf(stderr,"\nWARNING: something is wrong with the energy sums, will not use exact averages\n");
 832             ener_old->nsum_prev = 0;
 833         }
 834         else
 835         {
 836             ener_old->nsum_prev = nstep_all;
 837         }
 838         /* Copy all sums to ener_prev */
 839         for(i=0; i<fr->nre; i++)
 840         {
 841             ener_old->ener_prev[i].esum = fr->ener[i].esum;
 842             ener_old->ener_prev[i].eav  = fr->ener[i].eav;
 843         }
 844     }
 845
 846     ener_old->step_prev = fr->step;
 847 }
 848
 849 gmx_bool do_enx(ener_file_t ef,t_enxframe *fr)
 850 {
 851     int       file_version=-1;
 852     int       i,b;
 853     gmx_bool      bRead,bOK,bOK1,bSane;
 854     real      tmp1,tmp2,rdum;
 855     char      buf[22];
 856     /*int       d_size;*/
 857
 858     bOK = TRUE;
 859     bRead = gmx_fio_getread(ef->fio);
 860     if (!bRead)
 861     {
 862         fr->e_size = fr->nre*sizeof(fr->ener[0].e)*4;
 863         /*d_size = fr->ndisre*(sizeof(real)*2);*/
 864     }
 865     gmx_fio_checktype(ef->fio);
 866
 867     if (!do_eheader(ef,&file_version,fr,-1,NULL,&bOK))
 868     {
 869         if (bRead)
 870         {
 871             fprintf(stderr,"\rLast energy frame read %d time %8.3f         ",
 872                     ef->framenr-1,ef->frametime);
 873             if (!bOK)
 874             {
 875                 fprintf(stderr,
 876                         "\nWARNING: Incomplete energy frame: nr %d time %8.3f\n",
 877                         ef->framenr,fr->t);
 878             }
 879         }
 880         else
 881         {
 882             gmx_file("Cannot write energy file header; maybe you are out of disk space?");
 883         }
 884         return FALSE;
 885     }
 886     if (bRead)
 887     {
 888         if ((ef->framenr <   20 || ef->framenr %   10 == 0) &&
 889             (ef->framenr <  200 || ef->framenr %  100 == 0) &&
 890             (ef->framenr < 2000 || ef->framenr % 1000 == 0))
 891         {
 892             fprintf(stderr,"\rReading energy frame %6d time %8.3f         ",
 893                     ef->framenr,fr->t);
 894         }
 895         ef->framenr++;
 896         ef->frametime = fr->t;
 897     }
 898     /* Check sanity of this header */
 899     bSane = fr->nre > 0 ;
 900     for(b=0; b<fr->nblock; b++)
 901     {
 902         bSane = bSane || (fr->block[b].nsub > 0);
 903     }
 904     if (!((fr->step >= 0) && bSane))
 905     {
 906         fprintf(stderr,"\nWARNING: there may be something wrong with energy file %s\n",
 907                 gmx_fio_getname(ef->fio));
 908         fprintf(stderr,"Found: step=%s, nre=%d, nblock=%d, time=%g.\n"
 909                 "Trying to skip frame expect a crash though\n",
 910                 gmx_step_str(fr->step,buf),fr->nre,fr->nblock,fr->t);
 911     }
 912     if (bRead && fr->nre > fr->e_alloc)
 913     {
 914         srenew(fr->ener,fr->nre);
 915         for(i=fr->e_alloc; (i<fr->nre); i++)
 916         {
 917             fr->ener[i].e    = 0;
 918             fr->ener[i].eav  = 0;
 919             fr->ener[i].esum = 0;
 920         }
 921         fr->e_alloc = fr->nre;
 922     }
 923
 924     for(i=0; i<fr->nre; i++)
 925     {
 926         bOK = bOK && gmx_fio_do_real(ef->fio, fr->ener[i].e);
 927
 928         /* Do not store sums of length 1,
 929          * since this does not add information.
 930          */
 931         if (file_version == 1 ||
 932             (bRead && fr->nsum > 0) || fr->nsum > 1)
 933         {
 934             tmp1 = fr->ener[i].eav;
 935             bOK = bOK && gmx_fio_do_real(ef->fio, tmp1);
 936             if (bRead)
 937                 fr->ener[i].eav = tmp1;
 938
 939             /* This is to save only in single precision (unless compiled in DP) */
 940             tmp2 = fr->ener[i].esum;
 941             bOK = bOK && gmx_fio_do_real(ef->fio, tmp2);
 942             if (bRead)
 943                 fr->ener[i].esum = tmp2;
 944
 945             if (file_version == 1)
 946             {
 947                 /* Old, unused real */
 948                 rdum = 0;
 949                 bOK = bOK && gmx_fio_do_real(ef->fio, rdum);
 950             }
 951         }
 952     }
 953
 954     /* Here we can not check for file_version==1, since one could have
 955      * continued an old format simulation with a new one with mdrun -append.
 956      */
 957     if (bRead && ef->eo.bOldFileOpen)
 958     {
 959         /* Convert old full simulation sums to sums between energy frames */
 960         convert_full_sums(&(ef->eo),fr);
 961     }
 962     /* read the blocks */
 963     for(b=0; b<fr->nblock; b++)
 964     {
 965         /* now read the subblocks. */
 966         int nsub=fr->block[b].nsub; /* shortcut */
 967         int i;
 968
 969         for(i=0;i<nsub;i++)
 970         {
 971             t_enxsubblock *sub=&(fr->block[b].sub[i]); /* shortcut */
 972
 973             if (bRead)
 974             {
 975                 enxsubblock_alloc(sub);
 976             }
 977
 978             /* read/write data */
 979             bOK1=TRUE;
 980             switch (sub->type)
 981             {
 982                 case xdr_datatype_float:
 983                     bOK1=gmx_fio_ndo_float(ef->fio, sub->fval, sub->nr);
 984                     break;
 985                 case xdr_datatype_double:
 986                     bOK1=gmx_fio_ndo_double(ef->fio, sub->dval, sub->nr);
 987                     break;
 988                 case xdr_datatype_int:
 989                     bOK1=gmx_fio_ndo_int(ef->fio, sub->ival, sub->nr);
 990                     break;
 991                 case xdr_datatype_large_int:
 992                     bOK1=gmx_fio_ndo_gmx_large_int(ef->fio, sub->lval, sub->nr);
 993                     break;
 994                 case xdr_datatype_char:
 995                     bOK1=gmx_fio_ndo_uchar(ef->fio, sub->cval, sub->nr);
 996                     break;
 997                 case xdr_datatype_string:
 998                     bOK1=gmx_fio_ndo_string(ef->fio, sub->sval, sub->nr);
 999                     break;
1000                 default:
1001                     gmx_incons("Reading unknown block data type: this file is corrupted or from the future");
1002             }
1003             bOK = bOK && bOK1;
1004         }
1005     }
1006
1007     if(!bRead)
1008     {
1009         if( gmx_fio_flush(ef->fio) != 0)
1010         {
1011             gmx_file("Cannot write energy file; maybe you are out of disk space?");
1012         }
1013     }
1014
1015     if (!bOK)
1016     {
1017         if (bRead)
1018         {
1019             fprintf(stderr,"\nLast energy frame read %d",
1020                     ef->framenr-1);
1021             fprintf(stderr,"\nWARNING: Incomplete energy frame: nr %d time %8.3f\n",
1022                     ef->framenr,fr->t);
1023         }
1024         else
1025         {
1026             gmx_fatal(FARGS,"could not write energies");
1027         }
1028         return FALSE;
1029     }
1030
1031     return TRUE;
1032 }
1033
1034 static real find_energy(const char *name, int nre, gmx_enxnm_t *enm,
1035                         t_enxframe *fr)
1036 {
1037     int i;
1038
1039     for(i=0; i<nre; i++)
1040     {
1041         if (strcmp(enm[i].name,name) == 0)
1042         {
1043             return  fr->ener[i].e;
1044         }
1045     }
1046
1047     gmx_fatal(FARGS,"Could not find energy term named '%s'",name);
1048
1049     return 0;
1050 }
1051
1052
1053 void get_enx_state(const char *fn, real t, gmx_groups_t *groups, t_inputrec *ir,
1054                    t_state *state)
1055 {
1056   /* Should match the names in mdebin.c */
1057   static const char *boxvel_nm[] = {
1058   "Box-Vel-XX", "Box-Vel-YY", "Box-Vel-ZZ",
1059   "Box-Vel-YX", "Box-Vel-ZX", "Box-Vel-ZY"
1060   };
1061
1062   static const char *pcouplmu_nm[] = {
1063     "Pcoupl-Mu-XX", "Pcoupl-Mu-YY", "Pcoupl-Mu-ZZ",
1064     "Pcoupl-Mu-YX", "Pcoupl-Mu-ZX", "Pcoupl-Mu-ZY"
1065   };
1066   static const char *baro_nm[] = {
1067     "Barostat"
1068   };
1069
1070
1071   int ind0[] = { XX,YY,ZZ,YY,ZZ,ZZ };
1072   int ind1[] = { XX,YY,ZZ,XX,XX,YY };
1073   int nre,nfr,i,j,ni,npcoupl;
1074   char       buf[STRLEN];
1075   const char *bufi;
1076   gmx_enxnm_t *enm=NULL;
1077   t_enxframe *fr;
1078   ener_file_t in;
1079
1080   in = open_enx(fn,"r");
1081   do_enxnms(in,&nre,&enm);
1082   snew(fr,1);
1083   nfr = 0;
1084   while ((nfr==0 || fr->t != t) && do_enx(in,fr)) {
1085     nfr++;
1086   }
1087   close_enx(in);
1088   fprintf(stderr,"\n");
1089
1090   if (nfr == 0 || fr->t != t)
1091     gmx_fatal(FARGS,"Could not find frame with time %f in '%s'",t,fn);
1092
1093   npcoupl = TRICLINIC(ir->compress) ? 6 : 3;
1094   if (ir->epc == epcPARRINELLORAHMAN) {
1095     clear_mat(state->boxv);
1096     for(i=0; i<npcoupl; i++) {
1097       state->boxv[ind0[i]][ind1[i]] =
1098         find_energy(boxvel_nm[i],nre,enm,fr);
1099     }
1100     fprintf(stderr,"\nREAD %d BOX VELOCITIES FROM %s\n\n",npcoupl,fn);
1101   }
1102
1103   if (ir->etc == etcNOSEHOOVER)
1104   {
1105       char cns[20];
1106
1107       cns[0] = '\0';
1108
1109       for(i=0; i<state->ngtc; i++) {
1110           ni = groups->grps[egcTC].nm_ind[i];
1111           bufi = *(groups->grpname[ni]);
1112           for(j=0; (j<state->nhchainlength); j++)
1113           {
1114               if (IR_NVT_TROTTER(ir))
1115               {
1116                   sprintf(cns,"-%d",j);
1117               }
1118               sprintf(buf,"Xi%s-%s",cns,bufi);
1119               state->nosehoover_xi[i] = find_energy(buf,nre,enm,fr);
1120               sprintf(buf,"vXi%s-%s",cns,bufi);
1121               state->nosehoover_vxi[i] = find_energy(buf,nre,enm,fr);
1122           }
1123
1124       }
1125       fprintf(stderr,"\nREAD %d NOSE-HOOVER Xi chains FROM %s\n\n",state->ngtc,fn);
1126
1127       if (IR_NPT_TROTTER(ir) || IR_NPH_TROTTER(ir))
1128       {
1129           for(i=0; i<state->nnhpres; i++) {
1130               bufi = baro_nm[0]; /* All barostat DOF's together for now */
1131               for(j=0; (j<state->nhchainlength); j++)
1132               {
1133                   sprintf(buf,"Xi-%d-%s",j,bufi);
1134                   state->nhpres_xi[i] = find_energy(buf,nre,enm,fr);
1135                   sprintf(buf,"vXi-%d-%s",j,bufi);
1136                   state->nhpres_vxi[i] = find_energy(buf,nre,enm,fr);
1137               }
1138           }
1139           fprintf(stderr,"\nREAD %d NOSE-HOOVER BAROSTAT Xi chains FROM %s\n\n",state->nnhpres,fn);
1140       }
1141   }
1142
1143   free_enxnms(nre,enm);
1144   free_enxframe(fr);
1145   sfree(fr);
1146 }
1147