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 disk space?");
 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 disk space?");
 366             }
 367         }
 368         else
 369         {
 370             nm->unit = strdup("kJ/mol");
 371         }
 372     }
 373 }
 374
 375 void do_enxnms(ener_file_t ef,int *nre,gmx_enxnm_t **nms)
 376 {
 377     int  magic=-55555;
 378     XDR  *xdr;
 379     gmx_bool bRead = gmx_fio_getread(ef->fio);
 380     int  file_version;
 381     int  i;
 382
 383     gmx_fio_checktype(ef->fio);
 384
 385     xdr = gmx_fio_getxdr(ef->fio);
 386
 387     if (!xdr_int(xdr,&magic))
 388     {
 389         if(!bRead)
 390         {
 391             gmx_file("Cannot write energy names to file; maybe you are out of disk space?");
 392         }
 393         *nre=0;
 394         return;
 395     }
 396     if (magic > 0)
 397     {
 398         /* Assume this is an old edr format */
 399         file_version = 1;
 400         *nre = magic;
 401         ef->eo.bOldFileOpen = TRUE;
 402         ef->eo.bReadFirstStep = FALSE;
 403         srenew(ef->eo.ener_prev,*nre);
 404     }
 405     else
 406     {
 407         ef->eo.bOldFileOpen=FALSE;
 408
 409         if (magic != -55555)
 410         {
 411             gmx_fatal(FARGS,"Energy names magic number mismatch, this is not a GROMACS edr file");
 412         }
 413         file_version = enx_version;
 414         xdr_int(xdr,&file_version);
 415         if (file_version > enx_version)
 416         {
 417             gmx_fatal(FARGS,"reading tpx file (%s) version %d with version %d program",gmx_fio_getname(ef->fio),file_version,enx_version);
 418         }
 419         xdr_int(xdr,nre);
 420     }
 421     if (file_version != enx_version)
 422     {
 423         fprintf(stderr,"Note: enx file_version %d, software version %d\n",
 424                 file_version,enx_version);
 425     }
 426
 427     edr_strings(xdr,bRead,file_version,*nre,nms);
 428 }
 429
 430 static gmx_bool do_eheader(ener_file_t ef,int *file_version,t_enxframe *fr,
 431                        int nre_test,gmx_bool *bWrongPrecision,gmx_bool *bOK)
 432 {
 433     int  magic=-7777777;
 434     real first_real_to_check;
 435     int  b,i,zero=0,dum=0;
 436     gmx_bool bRead = gmx_fio_getread(ef->fio);
 437     int  tempfix_nr=0;
 438     int  ndisre=0;
 439     int  startb=0;
 440 #ifndef GMX_DOUBLE
 441     xdr_datatype dtreal=xdr_datatype_float;
 442 #else
 443     xdr_datatype dtreal=xdr_datatype_double;
 444 #endif
 445
 446     if (bWrongPrecision)
 447     {
 448         *bWrongPrecision = FALSE;
 449     }
 450
 451     *bOK=TRUE;
 452     /* The original energy frame started with a real,
 453      * so we have to use a real for compatibility.
 454      * This is VERY DIRTY code, since do_eheader can be called
 455      * with the wrong precision set and then we could read r > -1e10,
 456      * while actually the intention was r < -1e10.
 457      * When nre_test >= 0, do_eheader should therefore terminate
 458      * before the number of i/o calls starts depending on what has been read
 459      * (which is the case for for instance the block sizes for variable
 460      * number of blocks, where this number is read before).
 461      */
 462     first_real_to_check = -2e10;
 463     if (!gmx_fio_do_real(ef->fio, first_real_to_check))
 464     {
 465         return FALSE;
 466     }
 467     if (first_real_to_check > -1e10)
 468     {
 469         /* Assume we are reading an old format */
 470         *file_version = 1;
 471         fr->t = first_real_to_check;
 472         if (!gmx_fio_do_int(ef->fio, dum))   *bOK = FALSE;
 473         fr->step = dum;
 474     }
 475     else
 476     {
 477         if (!gmx_fio_do_int(ef->fio, magic))       *bOK = FALSE;
 478         if (magic != -7777777)
 479         {
 480             enx_warning("Energy header magic number mismatch, this is not a GROMACS edr file");
 481             *bOK=FALSE;
 482             return FALSE;
 483         }
 484         *file_version = enx_version;
 485         if (!gmx_fio_do_int(ef->fio, *file_version)) *bOK = FALSE;
 486         if (*bOK && *file_version > enx_version)
 487         {
 488             gmx_fatal(FARGS,"reading tpx file (%s) version %d with version %d program",gmx_fio_getname(ef->fio),file_version,enx_version);
 489         }
 490         if (!gmx_fio_do_double(ef->fio, fr->t))       *bOK = FALSE;
 491         if (!gmx_fio_do_gmx_large_int(ef->fio, fr->step)) *bOK = FALSE;
 492         if (!bRead && fr->nsum == 1) {
 493             /* Do not store sums of length 1,
 494              * since this does not add information.
 495              */
 496             if (!gmx_fio_do_int(ef->fio, zero))      *bOK = FALSE;
 497         } else {
 498             if (!gmx_fio_do_int(ef->fio, fr->nsum))  *bOK = FALSE;
 499         }
 500         if (*file_version >= 3)
 501         {
 502             if (!gmx_fio_do_gmx_large_int(ef->fio, fr->nsteps)) *bOK = FALSE;
 503         }
 504         else
 505         {
 506             fr->nsteps = max(1,fr->nsum);
 507         }
 508         if (*file_version >= 5)
 509         {
 510             if (!gmx_fio_do_double(ef->fio, fr->dt)) *bOK = FALSE;
 511         }
 512         else
 513         {
 514             fr->dt = 0;
 515         }
 516     }
 517     if (!gmx_fio_do_int(ef->fio, fr->nre))     *bOK = FALSE;
 518     if (*file_version < 4)
 519     {
 520         if (!gmx_fio_do_int(ef->fio, ndisre))  *bOK = FALSE;
 521     }
 522     else
 523     {
 524         /* now reserved for possible future use */
 525         if (!gmx_fio_do_int(ef->fio, dum))  *bOK = FALSE;
 526     }
 527
 528     if (!gmx_fio_do_int(ef->fio, fr->nblock))  *bOK = FALSE;
 529     if (fr->nblock < 0) *bOK=FALSE;
 530
 531     if (ndisre!=0)
 532     {
 533         if (*file_version >= 4)
 534         {
 535             enx_warning("Distance restraint blocks in old style in new style file");
 536             *bOK=FALSE;
 537             return FALSE;
 538         }
 539         fr->nblock+=1;
 540     }
 541
 542
 543     /* Frames could have nre=0, so we can not rely only on the fr->nre check */
 544     if (bRead && nre_test >= 0 &&
 545         ((fr->nre > 0 && fr->nre != nre_test) ||
 546          fr->nre < 0 || ndisre < 0 || fr->nblock < 0))
 547     {
 548         *bWrongPrecision = TRUE;
 549         return *bOK;
 550     }
 551
 552     /* we now know what these should be, or we've already bailed out because
 553        of wrong precision */
 554     if ( *file_version==1 && (fr->t < 0 || fr->t > 1e20 || fr->step < 0 ) )
 555     {
 556         enx_warning("edr file with negative step number or unreasonable time (and without version number).");
 557         *bOK=FALSE;
 558         return FALSE;
 559     }
 560
 561
 562     if (*bOK && bRead)
 563     {
 564         add_blocks_enxframe(fr, fr->nblock);
 565     }
 566
 567     startb=0;
 568     if (ndisre>0)
 569     {
 570         /* sub[0] is the instantaneous data, sub[1] is time averaged */
 571         add_subblocks_enxblock(&(fr->block[0]), 2);
 572         fr->block[0].id=enxDISRE;
 573         fr->block[0].sub[0].nr=ndisre;
 574         fr->block[0].sub[1].nr=ndisre;
 575         fr->block[0].sub[0].type=dtreal;
 576         fr->block[0].sub[1].type=dtreal;
 577         startb++;
 578     }
 579
 580     /* read block header info */
 581     for(b=startb; b<fr->nblock; b++)
 582     {
 583         if (*file_version<4)
 584         {
 585             /* blocks in old version files always have 1 subblock that
 586                consists of reals. */
 587             int nrint;
 588
 589             if (bRead)
 590             {
 591                 add_subblocks_enxblock(&(fr->block[b]), 1);
 592             }
 593             else
 594             {
 595                 if (fr->block[b].nsub != 1)
 596                 {
 597                     gmx_incons("Writing an old version .edr file with too many subblocks");
 598                 }
 599                 if (fr->block[b].sub[0].type != dtreal)
 600                 {
 601                     gmx_incons("Writing an old version .edr file the wrong subblock type");
 602                 }
 603             }
 604             nrint = fr->block[b].sub[0].nr;
 605
 606             if (!gmx_fio_do_int(ef->fio, nrint))
 607             {
 608                 *bOK = FALSE;
 609             }
 610             fr->block[b].id          = b - startb;
 611             fr->block[b].sub[0].nr   = nrint;
 612             fr->block[b].sub[0].type = dtreal;
 613         }
 614         else
 615         {
 616             int i;
 617             /* in the new version files, the block header only contains
 618                the ID and the number of subblocks */
 619             int nsub=fr->block[b].nsub;
 620             *bOK = *bOK && gmx_fio_do_int(ef->fio, fr->block[b].id);
 621             *bOK = *bOK && gmx_fio_do_int(ef->fio, nsub);
 622
 623             fr->block[b].nsub=nsub;
 624             if (bRead)
 625                 add_subblocks_enxblock(&(fr->block[b]), nsub);
 626
 627             /* read/write type & size for each subblock */
 628             for(i=0;i<nsub;i++)
 629             {
 630                 t_enxsubblock *sub=&(fr->block[b].sub[i]); /* shortcut */
 631                 int typenr=sub->type;
 632
 633                 *bOK=*bOK && gmx_fio_do_int(ef->fio, typenr);
 634                 *bOK=*bOK && gmx_fio_do_int(ef->fio, sub->nr);
 635
 636                 sub->type = (xdr_datatype)typenr;
 637             }
 638         }
 639     }
 640     if (!gmx_fio_do_int(ef->fio, fr->e_size))  *bOK = FALSE;
 641
 642     /* now reserved for possible future use */
 643     if (!gmx_fio_do_int(ef->fio, dum))  *bOK = FALSE;
 644
 645     /* Do a dummy int to keep the format compatible with the old code */
 646     if (!gmx_fio_do_int(ef->fio, dum))         *bOK = FALSE;
 647
 648     if (*bOK && *file_version == 1 && nre_test < 0)
 649     {
 650 #if 0
 651         if (fp >= ener_old_nalloc)
 652         {
 653             gmx_incons("Problem with reading old format energy files");
 654         }
 655 #endif
 656
 657         if (!ef->eo.bReadFirstStep)
 658         {
 659             ef->eo.bReadFirstStep = TRUE;
 660             ef->eo.first_step     = fr->step;
 661             ef->eo.step_prev      = fr->step;
 662             ef->eo.nsum_prev      = 0;
 663         }
 664
 665         fr->nsum   = fr->step - ef->eo.first_step + 1;
 666         fr->nsteps = fr->step - ef->eo.step_prev;
 667         fr->dt     = 0;
 668     }
 669
 670     return *bOK;
 671 }
 672
 673 void free_enxnms(int n,gmx_enxnm_t *nms)
 674 {
 675     int i;
 676
 677     for(i=0; i<n; i++)
 678     {
 679         sfree(nms[i].name);
 680         sfree(nms[i].unit);
 681     }
 682
 683     sfree(nms);
 684 }
 685
 686 void close_enx(ener_file_t ef)
 687 {
 688     if(gmx_fio_close(ef->fio) != 0)
 689     {
 690         gmx_file("Cannot close energy file; it might be corrupt, or maybe you are out of disk space?");
 691     }
 692 }
 693
 694 static gmx_bool empty_file(const char *fn)
 695 {
 696     FILE *fp;
 697     char dum;
 698     int  ret;
 699     gmx_bool bEmpty;
 700
 701     fp = gmx_fio_fopen(fn,"r");
 702     ret = fread(&dum,sizeof(dum),1,fp);
 703     bEmpty = feof(fp);
 704     gmx_fio_fclose(fp);
 705
 706     return bEmpty;
 707 }
 708
 709
 710 ener_file_t open_enx(const char *fn,const char *mode)
 711 {
 712     int        nre,i;
 713     gmx_enxnm_t *nms=NULL;
 714     int        file_version=-1;
 715     t_enxframe *fr;
 716     gmx_bool       bWrongPrecision,bOK=TRUE;
 717     struct ener_file *ef;
 718
 719     snew(ef,1);
 720
 721     if (mode[0]=='r') {
 722         ef->fio=gmx_fio_open(fn,mode);
 723         gmx_fio_checktype(ef->fio);
 724         gmx_fio_setprecision(ef->fio,FALSE);
 725         do_enxnms(ef,&nre,&nms);
 726         snew(fr,1);
 727         do_eheader(ef,&file_version,fr,nre,&bWrongPrecision,&bOK);
 728         if(!bOK)
 729         {
 730             gmx_file("Cannot read energy file header. Corrupt file?");
 731         }
 732
 733         /* Now check whether this file is in single precision */
 734         if (!bWrongPrecision &&
 735             ((fr->e_size && (fr->nre == nre) &&
 736               (nre*4*(long int)sizeof(float) == fr->e_size)) ) )
 737         {
 738             fprintf(stderr,"Opened %s as single precision energy file\n",fn);
 739             free_enxnms(nre,nms);
 740         }
 741         else
 742         {
 743             gmx_fio_rewind(ef->fio);
 744             gmx_fio_checktype(ef->fio);
 745             gmx_fio_setprecision(ef->fio,TRUE);
 746             do_enxnms(ef,&nre,&nms);
 747             do_eheader(ef,&file_version,fr,nre,&bWrongPrecision,&bOK);
 748             if(!bOK)
 749             {
 750                 gmx_file("Cannot write energy file header; maybe you are out of disk space?");
 751             }
 752
 753             if (((fr->e_size && (fr->nre == nre) &&
 754                             (nre*4*(long int)sizeof(double) == fr->e_size)) ))
 755                 fprintf(stderr,"Opened %s as double precision energy file\n",
 756                         fn);
 757             else {
 758                 if (empty_file(fn))
 759                     gmx_fatal(FARGS,"File %s is empty",fn);
 760                 else
 761                     gmx_fatal(FARGS,"Energy file %s not recognized, maybe different CPU?",
 762                               fn);
 763             }
 764             free_enxnms(nre,nms);
 765         }
 766         free_enxframe(fr);
 767         sfree(fr);
 768         gmx_fio_rewind(ef->fio);
 769     }
 770     else
 771         ef->fio = gmx_fio_open(fn,mode);
 772
 773     ef->framenr=0;
 774     ef->frametime=0;
 775     return ef;
 776 }
 777
 778 t_fileio *enx_file_pointer(const ener_file_t ef)
 779 {
 780     return ef->fio;
 781 }
 782
 783 static void convert_full_sums(ener_old_t *ener_old,t_enxframe *fr)
 784 {
 785     int nstep_all;
 786     int ne,ns,i;
 787     double esum_all,eav_all;
 788
 789     if (fr->nsum > 0)
 790     {
 791         ne = 0;
 792         ns = 0;
 793         for(i=0; i<fr->nre; i++)
 794         {
 795             if (fr->ener[i].e    != 0) ne++;
 796             if (fr->ener[i].esum != 0) ns++;
 797         }
 798         if (ne > 0 && ns == 0)
 799         {
 800             /* We do not have all energy sums */
 801             fr->nsum = 0;
 802         }
 803     }
 804
 805     /* Convert old full simulation sums to sums between energy frames */
 806     nstep_all = fr->step - ener_old->first_step + 1;
 807     if (fr->nsum > 1 && fr->nsum == nstep_all && ener_old->nsum_prev > 0)
 808     {
 809         /* Set the new sum length: the frame step difference */
 810         fr->nsum = fr->step - ener_old->step_prev;
 811         for(i=0; i<fr->nre; i++)
 812         {
 813             esum_all = fr->ener[i].esum;
 814             eav_all  = fr->ener[i].eav;
 815             fr->ener[i].esum = esum_all - ener_old->ener_prev[i].esum;
 816             fr->ener[i].eav  = eav_all  - ener_old->ener_prev[i].eav
 817                 - dsqr(ener_old->ener_prev[i].esum/(nstep_all - fr->nsum)
 818                        - esum_all/nstep_all)*
 819                 (nstep_all - fr->nsum)*nstep_all/(double)fr->nsum;
 820             ener_old->ener_prev[i].esum = esum_all;
 821             ener_old->ener_prev[i].eav  = eav_all;
 822         }
 823         ener_old->nsum_prev = nstep_all;
 824     }
 825     else if (fr->nsum > 0)
 826     {
 827         if (fr->nsum != nstep_all)
 828         {
 829             fprintf(stderr,"\nWARNING: something is wrong with the energy sums, will not use exact averages\n");
 830             ener_old->nsum_prev = 0;
 831         }
 832         else
 833         {
 834             ener_old->nsum_prev = nstep_all;
 835         }
 836         /* Copy all sums to ener_prev */
 837         for(i=0; i<fr->nre; i++)
 838         {
 839             ener_old->ener_prev[i].esum = fr->ener[i].esum;
 840             ener_old->ener_prev[i].eav  = fr->ener[i].eav;
 841         }
 842     }
 843
 844     ener_old->step_prev = fr->step;
 845 }
 846
 847 gmx_bool do_enx(ener_file_t ef,t_enxframe *fr)
 848 {
 849     int       file_version=-1;
 850     int       i,b;
 851     gmx_bool      bRead,bOK,bOK1,bSane;
 852     real      tmp1,tmp2,rdum;
 853     char      buf[22];
 854     /*int       d_size;*/
 855
 856     bOK = TRUE;
 857     bRead = gmx_fio_getread(ef->fio);
 858     if (!bRead)
 859     {
 860         fr->e_size = fr->nre*sizeof(fr->ener[0].e)*4;
 861         /*d_size = fr->ndisre*(sizeof(real)*2);*/
 862     }
 863     gmx_fio_checktype(ef->fio);
 864
 865     if (!do_eheader(ef,&file_version,fr,-1,NULL,&bOK))
 866     {
 867         if (bRead)
 868         {
 869             fprintf(stderr,"\rLast energy frame read %d time %8.3f         ",
 870                     ef->framenr-1,ef->frametime);
 871             if (!bOK)
 872             {
 873                 fprintf(stderr,
 874                         "\nWARNING: Incomplete energy frame: nr %d time %8.3f\n",
 875                         ef->framenr,fr->t);
 876             }
 877         }
 878         else
 879         {
 880             gmx_file("Cannot write energy file header; maybe you are out of disk space?");
 881         }
 882         return FALSE;
 883     }
 884     if (bRead)
 885     {
 886         if ((ef->framenr <   20 || ef->framenr %   10 == 0) &&
 887             (ef->framenr <  200 || ef->framenr %  100 == 0) &&
 888             (ef->framenr < 2000 || ef->framenr % 1000 == 0))
 889         {
 890             fprintf(stderr,"\rReading energy frame %6d time %8.3f         ",
 891                     ef->framenr,fr->t);
 892         }
 893         ef->framenr++;
 894         ef->frametime = fr->t;
 895     }
 896     /* Check sanity of this header */
 897     bSane = fr->nre > 0 ;
 898     for(b=0; b<fr->nblock; b++)
 899     {
 900         bSane = bSane || (fr->block[b].nsub > 0);
 901     }
 902     if (!((fr->step >= 0) && bSane))
 903     {
 904         fprintf(stderr,"\nWARNING: there may be something wrong with energy file %s\n",
 905                 gmx_fio_getname(ef->fio));
 906         fprintf(stderr,"Found: step=%s, nre=%d, nblock=%d, time=%g.\n"
 907                 "Trying to skip frame expect a crash though\n",
 908                 gmx_step_str(fr->step,buf),fr->nre,fr->nblock,fr->t);
 909     }
 910     if (bRead && fr->nre > fr->e_alloc)
 911     {
 912         srenew(fr->ener,fr->nre);
 913         for(i=fr->e_alloc; (i<fr->nre); i++)
 914         {
 915             fr->ener[i].e    = 0;
 916             fr->ener[i].eav  = 0;
 917             fr->ener[i].esum = 0;
 918         }
 919         fr->e_alloc = fr->nre;
 920     }
 921
 922     for(i=0; i<fr->nre; i++)
 923     {
 924         bOK = bOK && gmx_fio_do_real(ef->fio, fr->ener[i].e);
 925
 926         /* Do not store sums of length 1,
 927          * since this does not add information.
 928          */
 929         if (file_version == 1 ||
 930             (bRead && fr->nsum > 0) || fr->nsum > 1)
 931         {
 932             tmp1 = fr->ener[i].eav;
 933             bOK = bOK && gmx_fio_do_real(ef->fio, tmp1);
 934             if (bRead)
 935                 fr->ener[i].eav = tmp1;
 936
 937             /* This is to save only in single precision (unless compiled in DP) */
 938             tmp2 = fr->ener[i].esum;
 939             bOK = bOK && gmx_fio_do_real(ef->fio, tmp2);
 940             if (bRead)
 941                 fr->ener[i].esum = tmp2;
 942
 943             if (file_version == 1)
 944             {
 945                 /* Old, unused real */
 946                 rdum = 0;
 947                 bOK = bOK && gmx_fio_do_real(ef->fio, rdum);
 948             }
 949         }
 950     }
 951
 952     /* Here we can not check for file_version==1, since one could have
 953      * continued an old format simulation with a new one with mdrun -append.
 954      */
 955     if (bRead && ef->eo.bOldFileOpen)
 956     {
 957         /* Convert old full simulation sums to sums between energy frames */
 958         convert_full_sums(&(ef->eo),fr);
 959     }
 960     /* read the blocks */
 961     for(b=0; b<fr->nblock; b++)
 962     {
 963         /* now read the subblocks. */
 964         int nsub=fr->block[b].nsub; /* shortcut */
 965         int i;
 966
 967         for(i=0;i<nsub;i++)
 968         {
 969             t_enxsubblock *sub=&(fr->block[b].sub[i]); /* shortcut */
 970
 971             if (bRead)
 972             {
 973                 enxsubblock_alloc(sub);
 974             }
 975
 976             /* read/write data */
 977             bOK1=TRUE;
 978             switch (sub->type)
 979             {
 980                 case xdr_datatype_float:
 981                     bOK1=gmx_fio_ndo_float(ef->fio, sub->fval, sub->nr);
 982                     break;
 983                 case xdr_datatype_double:
 984                     bOK1=gmx_fio_ndo_double(ef->fio, sub->dval, sub->nr);
 985                     break;
 986                 case xdr_datatype_int:
 987                     bOK1=gmx_fio_ndo_int(ef->fio, sub->ival, sub->nr);
 988                     break;
 989                 case xdr_datatype_large_int:
 990                     bOK1=gmx_fio_ndo_gmx_large_int(ef->fio, sub->lval, sub->nr);
 991                     break;
 992                 case xdr_datatype_char:
 993                     bOK1=gmx_fio_ndo_uchar(ef->fio, sub->cval, sub->nr);
 994                     break;
 995                 case xdr_datatype_string:
 996                     bOK1=gmx_fio_ndo_string(ef->fio, sub->sval, sub->nr);
 997                     break;
 998                 default:
 999                     gmx_incons("Reading unknown block data type: this file is corrupted or from the future");
1000             }
1001             bOK = bOK && bOK1;
1002         }
1003     }
1004
1005     if(!bRead)
1006     {
1007         if( gmx_fio_flush(ef->fio) != 0)
1008         {
1009             gmx_file("Cannot write energy file; maybe you are out of disk space?");
1010         }
1011     }
1012
1013     if (!bOK)
1014     {
1015         if (bRead)
1016         {
1017             fprintf(stderr,"\nLast energy frame read %d",
1018                     ef->framenr-1);
1019             fprintf(stderr,"\nWARNING: Incomplete energy frame: nr %d time %8.3f\n",
1020                     ef->framenr,fr->t);
1021         }
1022         else
1023         {
1024             gmx_fatal(FARGS,"could not write energies");
1025         }
1026         return FALSE;
1027     }
1028
1029     return TRUE;
1030 }
1031
1032 static real find_energy(const char *name, int nre, gmx_enxnm_t *enm,
1033                         t_enxframe *fr)
1034 {
1035     int i;
1036
1037     for(i=0; i<nre; i++)
1038     {
1039         if (strcmp(enm[i].name,name) == 0)
1040         {
1041             return  fr->ener[i].e;
1042         }
1043     }
1044
1045     gmx_fatal(FARGS,"Could not find energy term named '%s'",name);
1046
1047     return 0;
1048 }
1049
1050
1051 void get_enx_state(const char *fn, real t, gmx_groups_t *groups, t_inputrec *ir,
1052                    t_state *state)
1053 {
1054   /* Should match the names in mdebin.c */
1055   static const char *boxvel_nm[] = {
1056   "Box-Vel-XX", "Box-Vel-YY", "Box-Vel-ZZ",
1057   "Box-Vel-YX", "Box-Vel-ZX", "Box-Vel-ZY"
1058   };
1059
1060   static const char *pcouplmu_nm[] = {
1061     "Pcoupl-Mu-XX", "Pcoupl-Mu-YY", "Pcoupl-Mu-ZZ",
1062     "Pcoupl-Mu-YX", "Pcoupl-Mu-ZX", "Pcoupl-Mu-ZY"
1063   };
1064   static const char *baro_nm[] = {
1065     "Barostat"
1066   };
1067
1068
1069   int ind0[] = { XX,YY,ZZ,YY,ZZ,ZZ };
1070   int ind1[] = { XX,YY,ZZ,XX,XX,YY };
1071   int nre,nfr,i,j,ni,npcoupl;
1072   char       buf[STRLEN];
1073   const char *bufi;
1074   gmx_enxnm_t *enm=NULL;
1075   t_enxframe *fr;
1076   ener_file_t in;
1077
1078   in = open_enx(fn,"r");
1079   do_enxnms(in,&nre,&enm);
1080   snew(fr,1);
1081   nfr = 0;
1082   while ((nfr==0 || fr->t != t) && do_enx(in,fr)) {
1083     nfr++;
1084   }
1085   close_enx(in);
1086   fprintf(stderr,"\n");
1087
1088   if (nfr == 0 || fr->t != t)
1089     gmx_fatal(FARGS,"Could not find frame with time %f in '%s'",t,fn);
1090
1091   npcoupl = TRICLINIC(ir->compress) ? 6 : 3;
1092   if (ir->epc == epcPARRINELLORAHMAN) {
1093     clear_mat(state->boxv);
1094     for(i=0; i<npcoupl; i++) {
1095       state->boxv[ind0[i]][ind1[i]] =
1096         find_energy(boxvel_nm[i],nre,enm,fr);
1097     }
1098     fprintf(stderr,"\nREAD %d BOX VELOCITIES FROM %s\n\n",npcoupl,fn);
1099   }
1100
1101   if (ir->etc == etcNOSEHOOVER)
1102   {
1103       char cns[20];
1104
1105       cns[0] = '\0';
1106
1107       for(i=0; i<state->ngtc; i++) {
1108           ni = groups->grps[egcTC].nm_ind[i];
1109           bufi = *(groups->grpname[ni]);
1110           for(j=0; (j<state->nhchainlength); j++)
1111           {
1112               if (IR_NVT_TROTTER(ir))
1113               {
1114                   sprintf(cns,"-%d",j);
1115               }
1116               sprintf(buf,"Xi%s-%s",cns,bufi);
1117               state->nosehoover_xi[i] = find_energy(buf,nre,enm,fr);
1118               sprintf(buf,"vXi%s-%s",cns,bufi);
1119               state->nosehoover_vxi[i] = find_energy(buf,nre,enm,fr);
1120           }
1121
1122       }
1123       fprintf(stderr,"\nREAD %d NOSE-HOOVER Xi chains FROM %s\n\n",state->ngtc,fn);
1124
1125       if (IR_NPT_TROTTER(ir) || IR_NPH_TROTTER(ir))
1126       {
1127           for(i=0; i<state->nnhpres; i++) {
1128               bufi = baro_nm[0]; /* All barostat DOF's together for now */
1129               for(j=0; (j<state->nhchainlength); j++)
1130               {
1131                   sprintf(buf,"Xi-%d-%s",j,bufi);
1132                   state->nhpres_xi[i] = find_energy(buf,nre,enm,fr);
1133                   sprintf(buf,"vXi-%d-%s",j,bufi);
1134                   state->nhpres_vxi[i] = find_energy(buf,nre,enm,fr);
1135               }
1136           }
1137           fprintf(stderr,"\nREAD %d NOSE-HOOVER BAROSTAT Xi chains FROM %s\n\n",state->nnhpres,fn);
1138       }
1139   }
1140
1141   free_enxnms(nre,enm);
1142   free_enxframe(fr);
1143   sfree(fr);
1144 }
1145