src/gmxlib/enxio.c

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