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,2015, 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 #include "gmxpre.h"
  38
  39 #include "enxio.h"
  40
  41 #include <stdlib.h>
  42 #include <string.h>
  43
  44 #include "gromacs/fileio/gmxfio.h"
  45 #include "gromacs/fileio/gmxfio-xdr.h"
  46 #include "gromacs/fileio/xdrf.h"
  47 #include "gromacs/legacyheaders/macros.h"
  48 #include "gromacs/math/vec.h"
  49 #include "gromacs/pbcutil/pbc.h"
  50 #include "gromacs/topology/topology.h"
  51 #include "gromacs/utility/fatalerror.h"
  52 #include "gromacs/utility/futil.h"
  53 #include "gromacs/utility/smalloc.h"
  54
  55 /* The source code in this file should be thread-safe.
  56          Please keep it that way. */
  57
  58 /* This number should be increased whenever the file format changes! */
  59 static const int enx_version = 5;
  60
  61 const char      *enx_block_id_name[] = {
  62     "Averaged orientation restraints",
  63     "Instantaneous orientation restraints",
  64     "Orientation restraint order tensor(s)",
  65     "Distance restraints",
  66     "Free energy data",
  67     "BAR histogram",
  68     "Delta H raw data"
  69 };
  70
  71
  72 /* Stuff for reading pre 4.1 energy files */
  73 typedef struct {
  74     gmx_bool     bOldFileOpen;   /* Is this an open old file? */
  75     gmx_bool     bReadFirstStep; /* Did we read the first step? */
  76     int          first_step;     /* First step in the energy file */
  77     int          step_prev;      /* Previous step */
  78     int          nsum_prev;      /* Previous step sum length */
  79     t_energy    *ener_prev;      /* Previous energy sums */
  80 } ener_old_t;
  81
  82 struct ener_file
  83 {
  84     ener_old_t eo;
  85     t_fileio  *fio;
  86     int        framenr;
  87     real       frametime;
  88 };
  89
  90 static void enxsubblock_init(t_enxsubblock *sb)
  91 {
  92     sb->nr = 0;
  93 #ifdef GMX_DOUBLE
  94     sb->type = xdr_datatype_double;
  95 #else
  96     sb->type = xdr_datatype_float;
  97 #endif
  98     sb->fval       = NULL;
  99     sb->dval       = NULL;
 100     sb->ival       = NULL;
 101     sb->lval       = NULL;
 102     sb->cval       = NULL;
 103     sb->sval       = NULL;
 104     sb->fval_alloc = 0;
 105     sb->dval_alloc = 0;
 106     sb->ival_alloc = 0;
 107     sb->lval_alloc = 0;
 108     sb->cval_alloc = 0;
 109     sb->sval_alloc = 0;
 110 }
 111
 112 static void enxsubblock_free(t_enxsubblock *sb)
 113 {
 114     if (sb->fval_alloc)
 115     {
 116         sfree(sb->fval);
 117         sb->fval_alloc = 0;
 118         sb->fval       = NULL;
 119     }
 120     if (sb->dval_alloc)
 121     {
 122         sfree(sb->dval);
 123         sb->dval_alloc = 0;
 124         sb->dval       = NULL;
 125     }
 126     if (sb->ival_alloc)
 127     {
 128         sfree(sb->ival);
 129         sb->ival_alloc = 0;
 130         sb->ival       = NULL;
 131     }
 132     if (sb->lval_alloc)
 133     {
 134         sfree(sb->lval);
 135         sb->lval_alloc = 0;
 136         sb->lval       = NULL;
 137     }
 138     if (sb->cval_alloc)
 139     {
 140         sfree(sb->cval);
 141         sb->cval_alloc = 0;
 142         sb->cval       = NULL;
 143     }
 144     if (sb->sval_alloc)
 145     {
 146         int i;
 147
 148         for (i = 0; i < sb->sval_alloc; i++)
 149         {
 150             if (sb->sval[i])
 151             {
 152                 sfree(sb->sval[i]);
 153             }
 154         }
 155         sfree(sb->sval);
 156         sb->sval_alloc = 0;
 157         sb->sval       = NULL;
 158     }
 159 }
 160
 161 /* allocate the appropriate amount of memory for the given type and nr */
 162 static void enxsubblock_alloc(t_enxsubblock *sb)
 163 {
 164     /* allocate the appropriate amount of memory */
 165     switch (sb->type)
 166     {
 167         case xdr_datatype_float:
 168             if (sb->nr > sb->fval_alloc)
 169             {
 170                 srenew(sb->fval, sb->nr);
 171                 sb->fval_alloc = sb->nr;
 172             }
 173             break;
 174         case xdr_datatype_double:
 175             if (sb->nr > sb->dval_alloc)
 176             {
 177                 srenew(sb->dval, sb->nr);
 178                 sb->dval_alloc = sb->nr;
 179             }
 180             break;
 181         case xdr_datatype_int:
 182             if (sb->nr > sb->ival_alloc)
 183             {
 184                 srenew(sb->ival, sb->nr);
 185                 sb->ival_alloc = sb->nr;
 186             }
 187             break;
 188         case xdr_datatype_int64:
 189             if (sb->nr > sb->lval_alloc)
 190             {
 191                 srenew(sb->lval, sb->nr);
 192                 sb->lval_alloc = sb->nr;
 193             }
 194             break;
 195         case xdr_datatype_char:
 196             if (sb->nr > sb->cval_alloc)
 197             {
 198                 srenew(sb->cval, sb->nr);
 199                 sb->cval_alloc = sb->nr;
 200             }
 201             break;
 202         case xdr_datatype_string:
 203             if (sb->nr > sb->sval_alloc)
 204             {
 205                 int i;
 206
 207                 srenew(sb->sval, sb->nr);
 208                 for (i = sb->sval_alloc; i < sb->nr; i++)
 209                 {
 210                     sb->sval[i] = NULL;
 211                 }
 212                 sb->sval_alloc = sb->nr;
 213             }
 214             break;
 215         default:
 216             gmx_incons("Unknown block type: this file is corrupted or from the future");
 217     }
 218 }
 219
 220 static void enxblock_init(t_enxblock *eb)
 221 {
 222     eb->id         = enxOR;
 223     eb->nsub       = 0;
 224     eb->sub        = NULL;
 225     eb->nsub_alloc = 0;
 226 }
 227
 228 static void enxblock_free(t_enxblock *eb)
 229 {
 230     if (eb->nsub_alloc > 0)
 231     {
 232         int i;
 233         for (i = 0; i < eb->nsub_alloc; i++)
 234         {
 235             enxsubblock_free(&(eb->sub[i]));
 236         }
 237         sfree(eb->sub);
 238         eb->nsub_alloc = 0;
 239         eb->sub        = NULL;
 240     }
 241 }
 242
 243 void init_enxframe(t_enxframe *fr)
 244 {
 245     fr->e_alloc = 0;
 246     fr->ener    = NULL;
 247
 248     /*fr->d_alloc=0;*/
 249     fr->ener = NULL;
 250
 251     /*fr->ndisre=0;*/
 252
 253     fr->nblock       = 0;
 254     fr->nblock_alloc = 0;
 255     fr->block        = NULL;
 256 }
 257
 258
 259 void free_enxframe(t_enxframe *fr)
 260 {
 261     int b;
 262
 263     if (fr->e_alloc)
 264     {
 265         sfree(fr->ener);
 266     }
 267     for (b = 0; b < fr->nblock_alloc; b++)
 268     {
 269         enxblock_free(&(fr->block[b]));
 270     }
 271     sfree(fr->block);
 272 }
 273
 274 void add_blocks_enxframe(t_enxframe *fr, int n)
 275 {
 276     fr->nblock = n;
 277     if (n > fr->nblock_alloc)
 278     {
 279         int b;
 280
 281         srenew(fr->block, n);
 282         for (b = fr->nblock_alloc; b < fr->nblock; b++)
 283         {
 284             enxblock_init(&(fr->block[b]));
 285         }
 286         fr->nblock_alloc = n;
 287     }
 288 }
 289
 290 t_enxblock *find_block_id_enxframe(t_enxframe *ef, int id, t_enxblock *prev)
 291 {
 292     gmx_off_t starti = 0;
 293     gmx_off_t i;
 294
 295     if (prev)
 296     {
 297         starti = (prev - ef->block) + 1;
 298     }
 299     for (i = starti; i < ef->nblock; i++)
 300     {
 301         if (ef->block[i].id == id)
 302         {
 303             return &(ef->block[i]);
 304         }
 305     }
 306     return NULL;
 307 }
 308
 309 void add_subblocks_enxblock(t_enxblock *eb, int n)
 310 {
 311     eb->nsub = n;
 312     if (eb->nsub > eb->nsub_alloc)
 313     {
 314         int b;
 315
 316         srenew(eb->sub, n);
 317         for (b = eb->nsub_alloc; b < n; b++)
 318         {
 319             enxsubblock_init(&(eb->sub[b]));
 320         }
 321         eb->nsub_alloc = n;
 322     }
 323 }
 324
 325 static void enx_warning(const char *msg)
 326 {
 327     if (getenv("GMX_ENX_NO_FATAL") != NULL)
 328     {
 329         gmx_warning(msg);
 330     }
 331     else
 332     {
 333         gmx_fatal(FARGS, "%s\n%s",
 334                   msg,
 335                   "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");
 336     }
 337 }
 338
 339 static void edr_strings(XDR *xdr, gmx_bool bRead, int file_version,
 340                         int n, gmx_enxnm_t **nms)
 341 {
 342     int          i;
 343     gmx_enxnm_t *nm;
 344
 345     if (*nms == NULL)
 346     {
 347         snew(*nms, n);
 348     }
 349     for (i = 0; i < n; i++)
 350     {
 351         nm = &(*nms)[i];
 352         if (bRead)
 353         {
 354             if (nm->name)
 355             {
 356                 sfree(nm->name);
 357                 nm->name = NULL;
 358             }
 359             if (nm->unit)
 360             {
 361                 sfree(nm->unit);
 362                 nm->unit = NULL;
 363             }
 364         }
 365         if (!xdr_string(xdr, &(nm->name), STRLEN))
 366         {
 367             gmx_file("Cannot write energy names to file; maybe you are out of disk space?");
 368         }
 369         if (file_version >= 2)
 370         {
 371             if (!xdr_string(xdr, &(nm->unit), STRLEN))
 372             {
 373                 gmx_file("Cannot write energy names to file; maybe you are out of disk space?");
 374             }
 375         }
 376         else
 377         {
 378             nm->unit = gmx_strdup("kJ/mol");
 379         }
 380     }
 381 }
 382
 383 void do_enxnms(ener_file_t ef, int *nre, gmx_enxnm_t **nms)
 384 {
 385     int      magic = -55555;
 386     XDR     *xdr;
 387     gmx_bool bRead = gmx_fio_getread(ef->fio);
 388     int      file_version;
 389     int      i;
 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_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_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     /*int       d_size;*/
 923
 924     bOK   = TRUE;
 925     bRead = gmx_fio_getread(ef->fio);
 926     if (!bRead)
 927     {
 928         fr->e_size = fr->nre*sizeof(fr->ener[0].e)*4;
 929         /*d_size = fr->ndisre*(sizeof(real)*2);*/
 930     }
 931
 932     if (!do_eheader(ef, &file_version, fr, -1, NULL, &bOK))
 933     {
 934         if (bRead)
 935         {
 936             fprintf(stderr, "\rLast energy frame read %d time %8.3f         ",
 937                     ef->framenr-1, ef->frametime);
 938             if (!bOK)
 939             {
 940                 fprintf(stderr,
 941                         "\nWARNING: Incomplete energy frame: nr %d time %8.3f\n",
 942                         ef->framenr, fr->t);
 943             }
 944         }
 945         else
 946         {
 947             gmx_file("Cannot write energy file header; maybe you are out of disk space?");
 948         }
 949         return FALSE;
 950     }
 951     if (bRead)
 952     {
 953         if ((ef->framenr <   20 || ef->framenr %   10 == 0) &&
 954             (ef->framenr <  200 || ef->framenr %  100 == 0) &&
 955             (ef->framenr < 2000 || ef->framenr % 1000 == 0))
 956         {
 957             fprintf(stderr, "\rReading energy frame %6d time %8.3f         ",
 958                     ef->framenr, fr->t);
 959         }
 960         ef->framenr++;
 961         ef->frametime = fr->t;
 962     }
 963     /* Check sanity of this header */
 964     bSane = fr->nre > 0;
 965     for (b = 0; b < fr->nblock; b++)
 966     {
 967         bSane = bSane || (fr->block[b].nsub > 0);
 968     }
 969     if (!((fr->step >= 0) && bSane))
 970     {
 971         fprintf(stderr, "\nWARNING: there may be something wrong with energy file %s\n",
 972                 gmx_fio_getname(ef->fio));
 973         fprintf(stderr, "Found: step=%"GMX_PRId64 ", nre=%d, nblock=%d, time=%g.\n"
 974                 "Trying to skip frame expect a crash though\n",
 975                 fr->step, fr->nre, fr->nblock, fr->t);
 976     }
 977     if (bRead && fr->nre > fr->e_alloc)
 978     {
 979         srenew(fr->ener, fr->nre);
 980         for (i = fr->e_alloc; (i < fr->nre); i++)
 981         {
 982             fr->ener[i].e    = 0;
 983             fr->ener[i].eav  = 0;
 984             fr->ener[i].esum = 0;
 985         }
 986         fr->e_alloc = fr->nre;
 987     }
 988
 989     for (i = 0; i < fr->nre; i++)
 990     {
 991         bOK = bOK && gmx_fio_do_real(ef->fio, fr->ener[i].e);
 992
 993         /* Do not store sums of length 1,
 994          * since this does not add information.
 995          */
 996         if (file_version == 1 ||
 997             (bRead && fr->nsum > 0) || fr->nsum > 1)
 998         {
 999             tmp1 = fr->ener[i].eav;
1000             bOK  = bOK && gmx_fio_do_real(ef->fio, tmp1);
1001             if (bRead)
1002             {
1003                 fr->ener[i].eav = tmp1;
1004             }
1005
1006             /* This is to save only in single precision (unless compiled in DP) */
1007             tmp2 = fr->ener[i].esum;
1008             bOK  = bOK && gmx_fio_do_real(ef->fio, tmp2);
1009             if (bRead)
1010             {
1011                 fr->ener[i].esum = tmp2;
1012             }
1013
1014             if (file_version == 1)
1015             {
1016                 /* Old, unused real */
1017                 rdum = 0;
1018                 bOK  = bOK && gmx_fio_do_real(ef->fio, rdum);
1019             }
1020         }
1021     }
1022
1023     /* Here we can not check for file_version==1, since one could have
1024      * continued an old format simulation with a new one with mdrun -append.
1025      */
1026     if (bRead && ef->eo.bOldFileOpen)
1027     {
1028         /* Convert old full simulation sums to sums between energy frames */
1029         convert_full_sums(&(ef->eo), fr);
1030     }
1031     /* read the blocks */
1032     for (b = 0; b < fr->nblock; b++)
1033     {
1034         /* now read the subblocks. */
1035         int nsub = fr->block[b].nsub; /* shortcut */
1036         int i;
1037
1038         for (i = 0; i < nsub; i++)
1039         {
1040             t_enxsubblock *sub = &(fr->block[b].sub[i]); /* shortcut */
1041
1042             if (bRead)
1043             {
1044                 enxsubblock_alloc(sub);
1045             }
1046
1047             /* read/write data */
1048             bOK1 = TRUE;
1049             switch (sub->type)
1050             {
1051                 case xdr_datatype_float:
1052                     bOK1 = gmx_fio_ndo_float(ef->fio, sub->fval, sub->nr);
1053                     break;
1054                 case xdr_datatype_double:
1055                     bOK1 = gmx_fio_ndo_double(ef->fio, sub->dval, sub->nr);
1056                     break;
1057                 case xdr_datatype_int:
1058                     bOK1 = gmx_fio_ndo_int(ef->fio, sub->ival, sub->nr);
1059                     break;
1060                 case xdr_datatype_int64:
1061                     bOK1 = gmx_fio_ndo_int64(ef->fio, sub->lval, sub->nr);
1062                     break;
1063                 case xdr_datatype_char:
1064                     bOK1 = gmx_fio_ndo_uchar(ef->fio, sub->cval, sub->nr);
1065                     break;
1066                 case xdr_datatype_string:
1067                     bOK1 = gmx_fio_ndo_string(ef->fio, sub->sval, sub->nr);
1068                     break;
1069                 default:
1070                     gmx_incons("Reading unknown block data type: this file is corrupted or from the future");
1071             }
1072             bOK = bOK && bOK1;
1073         }
1074     }
1075
1076     if (!bRead)
1077     {
1078         if (gmx_fio_flush(ef->fio) != 0)
1079         {
1080             gmx_file("Cannot write energy file; maybe you are out of disk space?");
1081         }
1082     }
1083
1084     if (!bOK)
1085     {
1086         if (bRead)
1087         {
1088             fprintf(stderr, "\nLast energy frame read %d",
1089                     ef->framenr-1);
1090             fprintf(stderr, "\nWARNING: Incomplete energy frame: nr %d time %8.3f\n",
1091                     ef->framenr, fr->t);
1092         }
1093         else
1094         {
1095             gmx_fatal(FARGS, "could not write energies");
1096         }
1097         return FALSE;
1098     }
1099
1100     return TRUE;
1101 }
1102
1103 static real find_energy(const char *name, int nre, gmx_enxnm_t *enm,
1104                         t_enxframe *fr)
1105 {
1106     int i;
1107
1108     for (i = 0; i < nre; i++)
1109     {
1110         if (strcmp(enm[i].name, name) == 0)
1111         {
1112             return fr->ener[i].e;
1113         }
1114     }
1115
1116     gmx_fatal(FARGS, "Could not find energy term named '%s'", name);
1117
1118     return 0;
1119 }
1120
1121
1122 void get_enx_state(const char *fn, real t, gmx_groups_t *groups, t_inputrec *ir,
1123                    t_state *state)
1124 {
1125     /* Should match the names in mdebin.c */
1126     static const char *boxvel_nm[] = {
1127         "Box-Vel-XX", "Box-Vel-YY", "Box-Vel-ZZ",
1128         "Box-Vel-YX", "Box-Vel-ZX", "Box-Vel-ZY"
1129     };
1130
1131     static const char *pcouplmu_nm[] = {
1132         "Pcoupl-Mu-XX", "Pcoupl-Mu-YY", "Pcoupl-Mu-ZZ",
1133         "Pcoupl-Mu-YX", "Pcoupl-Mu-ZX", "Pcoupl-Mu-ZY"
1134     };
1135     static const char *baro_nm[] = {
1136         "Barostat"
1137     };
1138
1139
1140     int          ind0[] = { XX, YY, ZZ, YY, ZZ, ZZ };
1141     int          ind1[] = { XX, YY, ZZ, XX, XX, YY };
1142     int          nre, nfr, i, j, ni, npcoupl;
1143     char         buf[STRLEN];
1144     const char  *bufi;
1145     gmx_enxnm_t *enm = NULL;
1146     t_enxframe  *fr;
1147     ener_file_t  in;
1148
1149     in = open_enx(fn, "r");
1150     do_enxnms(in, &nre, &enm);
1151     snew(fr, 1);
1152     nfr = 0;
1153     while ((nfr == 0 || fr->t != t) && do_enx(in, fr))
1154     {
1155         nfr++;
1156     }
1157     close_enx(in);
1158     fprintf(stderr, "\n");
1159
1160     if (nfr == 0 || fr->t != t)
1161     {
1162         gmx_fatal(FARGS, "Could not find frame with time %f in '%s'", t, fn);
1163     }
1164
1165     npcoupl = TRICLINIC(ir->compress) ? 6 : 3;
1166     if (ir->epc == epcPARRINELLORAHMAN)
1167     {
1168         clear_mat(state->boxv);
1169         for (i = 0; i < npcoupl; i++)
1170         {
1171             state->boxv[ind0[i]][ind1[i]] =
1172                 find_energy(boxvel_nm[i], nre, enm, fr);
1173         }
1174         fprintf(stderr, "\nREAD %d BOX VELOCITIES FROM %s\n\n", npcoupl, fn);
1175     }
1176
1177     if (ir->etc == etcNOSEHOOVER)
1178     {
1179         char cns[20];
1180
1181         cns[0] = '\0';
1182
1183         for (i = 0; i < state->ngtc; i++)
1184         {
1185             ni   = groups->grps[egcTC].nm_ind[i];
1186             bufi = *(groups->grpname[ni]);
1187             for (j = 0; (j < state->nhchainlength); j++)
1188             {
1189                 if (IR_NVT_TROTTER(ir))
1190                 {
1191                     sprintf(cns, "-%d", j);
1192                 }
1193                 sprintf(buf, "Xi%s-%s", cns, bufi);
1194                 state->nosehoover_xi[i] = find_energy(buf, nre, enm, fr);
1195                 sprintf(buf, "vXi%s-%s", cns, bufi);
1196                 state->nosehoover_vxi[i] = find_energy(buf, nre, enm, fr);
1197             }
1198
1199         }
1200         fprintf(stderr, "\nREAD %d NOSE-HOOVER Xi chains FROM %s\n\n", state->ngtc, fn);
1201
1202         if (IR_NPT_TROTTER(ir) || IR_NPH_TROTTER(ir))
1203         {
1204             for (i = 0; i < state->nnhpres; i++)
1205             {
1206                 bufi = baro_nm[0]; /* All barostat DOF's together for now */
1207                 for (j = 0; (j < state->nhchainlength); j++)
1208                 {
1209                     sprintf(buf, "Xi-%d-%s", j, bufi);
1210                     state->nhpres_xi[i] = find_energy(buf, nre, enm, fr);
1211                     sprintf(buf, "vXi-%d-%s", j, bufi);
1212                     state->nhpres_vxi[i] = find_energy(buf, nre, enm, fr);
1213                 }
1214             }
1215             fprintf(stderr, "\nREAD %d NOSE-HOOVER BAROSTAT Xi chains FROM %s\n\n", state->nnhpres, fn);
1216         }
1217     }
1218
1219     free_enxnms(nre, enm);
1220     free_enxframe(fr);
1221     sfree(fr);
1222 }