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