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