src/gromacs/fileio/enxio.c

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