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