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 #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/xdrf.h"
  46 #include "gromacs/legacyheaders/macros.h"
  47 #include "gromacs/math/vec.h"
  48 #include "gromacs/pbcutil/pbc.h"
  49 #include "gromacs/topology/topology.h"
  50 #include "gromacs/utility/fatalerror.h"
  51 #include "gromacs/utility/futil.h"
  52 #include "gromacs/utility/smalloc.h"
  53
  54 /* The source code in this file should be thread-safe.
  55          Please keep it that way. */
  56
  57 /* This number should be increased whenever the file format changes! */
  58 static const int enx_version = 5;
  59
  60 const char      *enx_block_id_name[] = {
  61     "Averaged orientation restraints",
  62     "Instantaneous orientation restraints",
  63     "Orientation restraint order tensor(s)",
  64     "Distance restraints",
  65     "Free energy data",
  66     "BAR histogram",
  67     "Delta H raw data"
  68 };
  69
  70
  71 /* Stuff for reading pre 4.1 energy files */
  72 typedef struct {
  73     gmx_bool     bOldFileOpen;   /* Is this an open old file? */
  74     gmx_bool     bReadFirstStep; /* Did we read the first step? */
  75     int          first_step;     /* First step in the energy file */
  76     int          step_prev;      /* Previous step */
  77     int          nsum_prev;      /* Previous step sum length */
  78     t_energy    *ener_prev;      /* Previous energy sums */
  79 } ener_old_t;
  80
  81 struct ener_file
  82 {
  83     ener_old_t eo;
  84     t_fileio  *fio;
  85     int        framenr;
  86     real       frametime;
  87 };
  88
  89 static void enxsubblock_init(t_enxsubblock *sb)
  90 {
  91     sb->nr = 0;
  92 #ifdef GMX_DOUBLE
  93     sb->type = xdr_datatype_double;
  94 #else
  95     sb->type = xdr_datatype_float;
  96 #endif
  97     sb->fval       = NULL;
  98     sb->dval       = NULL;
  99     sb->ival       = NULL;
 100     sb->lval       = NULL;
 101     sb->cval       = NULL;
 102     sb->sval       = NULL;
 103     sb->fval_alloc = 0;
 104     sb->dval_alloc = 0;
 105     sb->ival_alloc = 0;
 106     sb->lval_alloc = 0;
 107     sb->cval_alloc = 0;
 108     sb->sval_alloc = 0;
 109 }
 110
 111 static void enxsubblock_free(t_enxsubblock *sb)
 112 {
 113     if (sb->fval_alloc)
 114     {
 115         sfree(sb->fval);
 116         sb->fval_alloc = 0;
 117         sb->fval       = NULL;
 118     }
 119     if (sb->dval_alloc)
 120     {
 121         sfree(sb->dval);
 122         sb->dval_alloc = 0;
 123         sb->dval       = NULL;
 124     }
 125     if (sb->ival_alloc)
 126     {
 127         sfree(sb->ival);
 128         sb->ival_alloc = 0;
 129         sb->ival       = NULL;
 130     }
 131     if (sb->lval_alloc)
 132     {
 133         sfree(sb->lval);
 134         sb->lval_alloc = 0;
 135         sb->lval       = NULL;
 136     }
 137     if (sb->cval_alloc)
 138     {
 139         sfree(sb->cval);
 140         sb->cval_alloc = 0;
 141         sb->cval       = NULL;
 142     }
 143     if (sb->sval_alloc)
 144     {
 145         int i;
 146
 147         for (i = 0; i < sb->sval_alloc; i++)
 148         {
 149             if (sb->sval[i])
 150             {
 151                 sfree(sb->sval[i]);
 152             }
 153         }
 154         sfree(sb->sval);
 155         sb->sval_alloc = 0;
 156         sb->sval       = NULL;
 157     }
 158 }
 159
 160 /* allocate the appropriate amount of memory for the given type and nr */
 161 static void enxsubblock_alloc(t_enxsubblock *sb)
 162 {
 163     /* allocate the appropriate amount of memory */
 164     switch (sb->type)
 165     {
 166         case xdr_datatype_float:
 167             if (sb->nr > sb->fval_alloc)
 168             {
 169                 srenew(sb->fval, sb->nr);
 170                 sb->fval_alloc = sb->nr;
 171             }
 172             break;
 173         case xdr_datatype_double:
 174             if (sb->nr > sb->dval_alloc)
 175             {
 176                 srenew(sb->dval, sb->nr);
 177                 sb->dval_alloc = sb->nr;
 178             }
 179             break;
 180         case xdr_datatype_int:
 181             if (sb->nr > sb->ival_alloc)
 182             {
 183                 srenew(sb->ival, sb->nr);
 184                 sb->ival_alloc = sb->nr;
 185             }
 186             break;
 187         case xdr_datatype_int64:
 188             if (sb->nr > sb->lval_alloc)
 189             {
 190                 srenew(sb->lval, sb->nr);
 191                 sb->lval_alloc = sb->nr;
 192             }
 193             break;
 194         case xdr_datatype_char:
 195             if (sb->nr > sb->cval_alloc)
 196             {
 197                 srenew(sb->cval, sb->nr);
 198                 sb->cval_alloc = sb->nr;
 199             }
 200             break;
 201         case xdr_datatype_string:
 202             if (sb->nr > sb->sval_alloc)
 203             {
 204                 int i;
 205
 206                 srenew(sb->sval, sb->nr);
 207                 for (i = sb->sval_alloc; i < sb->nr; i++)
 208                 {
 209                     sb->sval[i] = NULL;
 210                 }
 211                 sb->sval_alloc = sb->nr;
 212             }
 213             break;
 214         default:
 215             gmx_incons("Unknown block type: this file is corrupted or from the future");
 216     }
 217 }
 218
 219 static void enxblock_init(t_enxblock *eb)
 220 {
 221     eb->id         = enxOR;
 222     eb->nsub       = 0;
 223     eb->sub        = NULL;
 224     eb->nsub_alloc = 0;
 225 }
 226
 227 static void enxblock_free(t_enxblock *eb)
 228 {
 229     if (eb->nsub_alloc > 0)
 230     {
 231         int i;
 232         for (i = 0; i < eb->nsub_alloc; i++)
 233         {
 234             enxsubblock_free(&(eb->sub[i]));
 235         }
 236         sfree(eb->sub);
 237         eb->nsub_alloc = 0;
 238         eb->sub        = NULL;
 239     }
 240 }
 241
 242 void init_enxframe(t_enxframe *fr)
 243 {
 244     fr->e_alloc = 0;
 245     fr->ener    = NULL;
 246
 247     /*fr->d_alloc=0;*/
 248     fr->ener = NULL;
 249
 250     /*fr->ndisre=0;*/
 251
 252     fr->nblock       = 0;
 253     fr->nblock_alloc = 0;
 254     fr->block        = NULL;
 255 }
 256
 257
 258 void free_enxframe(t_enxframe *fr)
 259 {
 260     int b;
 261
 262     if (fr->e_alloc)
 263     {
 264         sfree(fr->ener);
 265     }
 266     for (b = 0; b < fr->nblock_alloc; b++)
 267     {
 268         enxblock_free(&(fr->block[b]));
 269     }
 270     sfree(fr->block);
 271 }
 272
 273 void add_blocks_enxframe(t_enxframe *fr, int n)
 274 {
 275     fr->nblock = n;
 276     if (n > fr->nblock_alloc)
 277     {
 278         int b;
 279
 280         srenew(fr->block, n);
 281         for (b = fr->nblock_alloc; b < fr->nblock; b++)
 282         {
 283             enxblock_init(&(fr->block[b]));
 284         }
 285         fr->nblock_alloc = n;
 286     }
 287 }
 288
 289 t_enxblock *find_block_id_enxframe(t_enxframe *ef, int id, t_enxblock *prev)
 290 {
 291     gmx_off_t starti = 0;
 292     gmx_off_t i;
 293
 294     if (prev)
 295     {
 296         starti = (prev - ef->block) + 1;
 297     }
 298     for (i = starti; i < ef->nblock; i++)
 299     {
 300         if (ef->block[i].id == id)
 301         {
 302             return &(ef->block[i]);
 303         }
 304     }
 305     return NULL;
 306 }
 307
 308 void add_subblocks_enxblock(t_enxblock *eb, int n)
 309 {
 310     eb->nsub = n;
 311     if (eb->nsub > eb->nsub_alloc)
 312     {
 313         int b;
 314
 315         srenew(eb->sub, n);
 316         for (b = eb->nsub_alloc; b < n; b++)
 317         {
 318             enxsubblock_init(&(eb->sub[b]));
 319         }
 320         eb->nsub_alloc = n;
 321     }
 322 }
 323
 324 static void enx_warning(const char *msg)
 325 {
 326     if (getenv("GMX_ENX_NO_FATAL") != NULL)
 327     {
 328         gmx_warning(msg);
 329     }
 330     else
 331     {
 332         gmx_fatal(FARGS, "%s\n%s",
 333                   msg,
 334                   "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");
 335     }
 336 }
 337
 338 static void edr_strings(XDR *xdr, gmx_bool bRead, int file_version,
 339                         int n, gmx_enxnm_t **nms)
 340 {
 341     int          i;
 342     gmx_enxnm_t *nm;
 343
 344     if (*nms == NULL)
 345     {
 346         snew(*nms, n);
 347     }
 348     for (i = 0; i < n; i++)
 349     {
 350         nm = &(*nms)[i];
 351         if (bRead)
 352         {
 353             if (nm->name)
 354             {
 355                 sfree(nm->name);
 356                 nm->name = NULL;
 357             }
 358             if (nm->unit)
 359             {
 360                 sfree(nm->unit);
 361                 nm->unit = NULL;
 362             }
 363         }
 364         if (!xdr_string(xdr, &(nm->name), STRLEN))
 365         {
 366             gmx_file("Cannot write energy names to file; maybe you are out of disk space?");
 367         }
 368         if (file_version >= 2)
 369         {
 370             if (!xdr_string(xdr, &(nm->unit), STRLEN))
 371             {
 372                 gmx_file("Cannot write energy names to file; maybe you are out of disk space?");
 373             }
 374         }
 375         else
 376         {
 377             nm->unit = gmx_strdup("kJ/mol");
 378         }
 379     }
 380 }
 381
 382 void do_enxnms(ener_file_t ef, int *nre, gmx_enxnm_t **nms)
 383 {
 384     int      magic = -55555;
 385     XDR     *xdr;
 386     gmx_bool bRead = gmx_fio_getread(ef->fio);
 387     int      file_version;
 388     int      i;
 389
 390     gmx_fio_checktype(ef->fio);
 391
 392     xdr = gmx_fio_getxdr(ef->fio);
 393
 394     if (!xdr_int(xdr, &magic))
 395     {
 396         if (!bRead)
 397         {
 398             gmx_file("Cannot write energy names to file; maybe you are out of disk space?");
 399         }
 400         *nre = 0;
 401         return;
 402     }
 403     if (magic > 0)
 404     {
 405         /* Assume this is an old edr format */
 406         file_version          = 1;
 407         *nre                  = magic;
 408         ef->eo.bOldFileOpen   = TRUE;
 409         ef->eo.bReadFirstStep = FALSE;
 410         srenew(ef->eo.ener_prev, *nre);
 411     }
 412     else
 413     {
 414         ef->eo.bOldFileOpen = FALSE;
 415
 416         if (magic != -55555)
 417         {
 418             gmx_fatal(FARGS, "Energy names magic number mismatch, this is not a GROMACS edr file");
 419         }
 420         file_version = enx_version;
 421         xdr_int(xdr, &file_version);
 422         if (file_version > enx_version)
 423         {
 424             gmx_fatal(FARGS, "reading tpx file (%s) version %d with version %d program", gmx_fio_getname(ef->fio), file_version, enx_version);
 425         }
 426         xdr_int(xdr, nre);
 427     }
 428     if (file_version != enx_version)
 429     {
 430         fprintf(stderr, "Note: enx file_version %d, software version %d\n",
 431                 file_version, enx_version);
 432     }
 433
 434     edr_strings(xdr, bRead, file_version, *nre, nms);
 435 }
 436
 437 static gmx_bool do_eheader(ener_file_t ef, int *file_version, t_enxframe *fr,
 438                            int nre_test, gmx_bool *bWrongPrecision, gmx_bool *bOK)
 439 {
 440     int          magic = -7777777;
 441     real         first_real_to_check;
 442     int          b, i, zero = 0, dum = 0;
 443     gmx_bool     bRead      = gmx_fio_getread(ef->fio);
 444     int          tempfix_nr = 0;
 445     int          ndisre     = 0;
 446     int          startb     = 0;
 447 #ifndef GMX_DOUBLE
 448     xdr_datatype dtreal = xdr_datatype_float;
 449 #else
 450     xdr_datatype dtreal = xdr_datatype_double;
 451 #endif
 452
 453     if (bWrongPrecision)
 454     {
 455         *bWrongPrecision = FALSE;
 456     }
 457
 458     *bOK = TRUE;
 459     /* The original energy frame started with a real,
 460      * so we have to use a real for compatibility.
 461      * This is VERY DIRTY code, since do_eheader can be called
 462      * with the wrong precision set and then we could read r > -1e10,
 463      * while actually the intention was r < -1e10.
 464      * When nre_test >= 0, do_eheader should therefore terminate
 465      * before the number of i/o calls starts depending on what has been read
 466      * (which is the case for for instance the block sizes for variable
 467      * number of blocks, where this number is read before).
 468      */
 469     first_real_to_check = -2e10;
 470     if (!gmx_fio_do_real(ef->fio, first_real_to_check))
 471     {
 472         return FALSE;
 473     }
 474     if (first_real_to_check > -1e10)
 475     {
 476         /* Assume we are reading an old format */
 477         *file_version = 1;
 478         fr->t         = first_real_to_check;
 479         if (!gmx_fio_do_int(ef->fio, dum))
 480         {
 481             *bOK = FALSE;
 482         }
 483         fr->step = dum;
 484     }
 485     else
 486     {
 487         if (!gmx_fio_do_int(ef->fio, magic))
 488         {
 489             *bOK = FALSE;
 490         }
 491         if (magic != -7777777)
 492         {
 493             enx_warning("Energy header magic number mismatch, this is not a GROMACS edr file");
 494             *bOK = FALSE;
 495             return FALSE;
 496         }
 497         *file_version = enx_version;
 498         if (!gmx_fio_do_int(ef->fio, *file_version))
 499         {
 500             *bOK = FALSE;
 501         }
 502         if (*bOK && *file_version > enx_version)
 503         {
 504             gmx_fatal(FARGS, "reading tpx file (%s) version %d with version %d program", gmx_fio_getname(ef->fio), file_version, enx_version);
 505         }
 506         if (!gmx_fio_do_double(ef->fio, fr->t))
 507         {
 508             *bOK = FALSE;
 509         }
 510         if (!gmx_fio_do_int64(ef->fio, fr->step))
 511         {
 512             *bOK = FALSE;
 513         }
 514         if (!bRead && fr->nsum == 1)
 515         {
 516             /* Do not store sums of length 1,
 517              * since this does not add information.
 518              */
 519             if (!gmx_fio_do_int(ef->fio, zero))
 520             {
 521                 *bOK = FALSE;
 522             }
 523         }
 524         else
 525         {
 526             if (!gmx_fio_do_int(ef->fio, fr->nsum))
 527             {
 528                 *bOK = FALSE;
 529             }
 530         }
 531         if (*file_version >= 3)
 532         {
 533             if (!gmx_fio_do_int64(ef->fio, fr->nsteps))
 534             {
 535                 *bOK = FALSE;
 536             }
 537         }
 538         else
 539         {
 540             fr->nsteps = max(1, fr->nsum);
 541         }
 542         if (*file_version >= 5)
 543         {
 544             if (!gmx_fio_do_double(ef->fio, fr->dt))
 545             {
 546                 *bOK = FALSE;
 547             }
 548         }
 549         else
 550         {
 551             fr->dt = 0;
 552         }
 553     }
 554     if (!gmx_fio_do_int(ef->fio, fr->nre))
 555     {
 556         *bOK = FALSE;
 557     }
 558     if (*file_version < 4)
 559     {
 560         if (!gmx_fio_do_int(ef->fio, ndisre))
 561         {
 562             *bOK = FALSE;
 563         }
 564     }
 565     else
 566     {
 567         /* now reserved for possible future use */
 568         if (!gmx_fio_do_int(ef->fio, dum))
 569         {
 570             *bOK = FALSE;
 571         }
 572     }
 573
 574     if (!gmx_fio_do_int(ef->fio, fr->nblock))
 575     {
 576         *bOK = FALSE;
 577     }
 578     if (fr->nblock < 0)
 579     {
 580         *bOK = FALSE;
 581     }
 582
 583     if (ndisre != 0)
 584     {
 585         if (*file_version >= 4)
 586         {
 587             enx_warning("Distance restraint blocks in old style in new style file");
 588             *bOK = FALSE;
 589             return FALSE;
 590         }
 591         fr->nblock += 1;
 592     }
 593
 594
 595     /* Frames could have nre=0, so we can not rely only on the fr->nre check */
 596     if (bRead && nre_test >= 0 &&
 597         ((fr->nre > 0 && fr->nre != nre_test) ||
 598          fr->nre < 0 || ndisre < 0 || fr->nblock < 0))
 599     {
 600         *bWrongPrecision = TRUE;
 601         return *bOK;
 602     }
 603
 604     /* we now know what these should be, or we've already bailed out because
 605        of wrong precision */
 606     if (*file_version == 1 && (fr->t < 0 || fr->t > 1e20 || fr->step < 0 ) )
 607     {
 608         enx_warning("edr file with negative step number or unreasonable time (and without version number).");
 609         *bOK = FALSE;
 610         return FALSE;
 611     }
 612
 613
 614     if (*bOK && bRead)
 615     {
 616         add_blocks_enxframe(fr, fr->nblock);
 617     }
 618
 619     startb = 0;
 620     if (ndisre > 0)
 621     {
 622         /* sub[0] is the instantaneous data, sub[1] is time averaged */
 623         add_subblocks_enxblock(&(fr->block[0]), 2);
 624         fr->block[0].id          = enxDISRE;
 625         fr->block[0].sub[0].nr   = ndisre;
 626         fr->block[0].sub[1].nr   = ndisre;
 627         fr->block[0].sub[0].type = dtreal;
 628         fr->block[0].sub[1].type = dtreal;
 629         startb++;
 630     }
 631
 632     /* read block header info */
 633     for (b = startb; b < fr->nblock; b++)
 634     {
 635         if (*file_version < 4)
 636         {
 637             /* blocks in old version files always have 1 subblock that
 638                consists of reals. */
 639             int nrint;
 640
 641             if (bRead)
 642             {
 643                 add_subblocks_enxblock(&(fr->block[b]), 1);
 644             }
 645             else
 646             {
 647                 if (fr->block[b].nsub != 1)
 648                 {
 649                     gmx_incons("Writing an old version .edr file with too many subblocks");
 650                 }
 651                 if (fr->block[b].sub[0].type != dtreal)
 652                 {
 653                     gmx_incons("Writing an old version .edr file the wrong subblock type");
 654                 }
 655             }
 656             nrint = fr->block[b].sub[0].nr;
 657
 658             if (!gmx_fio_do_int(ef->fio, nrint))
 659             {
 660                 *bOK = FALSE;
 661             }
 662             fr->block[b].id          = b - startb;
 663             fr->block[b].sub[0].nr   = nrint;
 664             fr->block[b].sub[0].type = dtreal;
 665         }
 666         else
 667         {
 668             int i;
 669             /* in the new version files, the block header only contains
 670                the ID and the number of subblocks */
 671             int nsub = fr->block[b].nsub;
 672             *bOK = *bOK && gmx_fio_do_int(ef->fio, fr->block[b].id);
 673             *bOK = *bOK && gmx_fio_do_int(ef->fio, nsub);
 674
 675             fr->block[b].nsub = nsub;
 676             if (bRead)
 677             {
 678                 add_subblocks_enxblock(&(fr->block[b]), nsub);
 679             }
 680
 681             /* read/write type & size for each subblock */
 682             for (i = 0; i < nsub; i++)
 683             {
 684                 t_enxsubblock *sub    = &(fr->block[b].sub[i]); /* shortcut */
 685                 int            typenr = sub->type;
 686
 687                 *bOK = *bOK && gmx_fio_do_int(ef->fio, typenr);
 688                 *bOK = *bOK && gmx_fio_do_int(ef->fio, sub->nr);
 689
 690                 sub->type = (xdr_datatype)typenr;
 691             }
 692         }
 693     }
 694     if (!gmx_fio_do_int(ef->fio, fr->e_size))
 695     {
 696         *bOK = FALSE;
 697     }
 698
 699     /* now reserved for possible future use */
 700     if (!gmx_fio_do_int(ef->fio, dum))
 701     {
 702         *bOK = FALSE;
 703     }
 704
 705     /* Do a dummy int to keep the format compatible with the old code */
 706     if (!gmx_fio_do_int(ef->fio, dum))
 707     {
 708         *bOK = FALSE;
 709     }
 710
 711     if (*bOK && *file_version == 1 && nre_test < 0)
 712     {
 713         if (!ef->eo.bReadFirstStep)
 714         {
 715             ef->eo.bReadFirstStep = TRUE;
 716             ef->eo.first_step     = fr->step;
 717             ef->eo.step_prev      = fr->step;
 718             ef->eo.nsum_prev      = 0;
 719         }
 720
 721         fr->nsum   = fr->step - ef->eo.first_step + 1;
 722         fr->nsteps = fr->step - ef->eo.step_prev;
 723         fr->dt     = 0;
 724     }
 725
 726     return *bOK;
 727 }
 728
 729 void free_enxnms(int n, gmx_enxnm_t *nms)
 730 {
 731     int i;
 732
 733     for (i = 0; i < n; i++)
 734     {
 735         sfree(nms[i].name);
 736         sfree(nms[i].unit);
 737     }
 738
 739     sfree(nms);
 740 }
 741
 742 void close_enx(ener_file_t ef)
 743 {
 744     if (gmx_fio_close(ef->fio) != 0)
 745     {
 746         gmx_file("Cannot close energy file; it might be corrupt, or maybe you are out of disk space?");
 747     }
 748 }
 749
 750 static gmx_bool empty_file(const char *fn)
 751 {
 752     FILE    *fp;
 753     char     dum;
 754     int      ret;
 755     gmx_bool bEmpty;
 756
 757     fp     = gmx_fio_fopen(fn, "r");
 758     ret    = fread(&dum, sizeof(dum), 1, fp);
 759     bEmpty = feof(fp);
 760     gmx_fio_fclose(fp);
 761
 762     return bEmpty;
 763 }
 764
 765
 766 ener_file_t open_enx(const char *fn, const char *mode)
 767 {
 768     int               nre, i;
 769     gmx_enxnm_t      *nms          = NULL;
 770     int               file_version = -1;
 771     t_enxframe       *fr;
 772     gmx_bool          bWrongPrecision, bOK = TRUE;
 773     struct ener_file *ef;
 774
 775     snew(ef, 1);
 776
 777     if (mode[0] == 'r')
 778     {
 779         ef->fio = gmx_fio_open(fn, mode);
 780         gmx_fio_checktype(ef->fio);
 781         gmx_fio_setprecision(ef->fio, FALSE);
 782         do_enxnms(ef, &nre, &nms);
 783         snew(fr, 1);
 784         do_eheader(ef, &file_version, fr, nre, &bWrongPrecision, &bOK);
 785         if (!bOK)
 786         {
 787             gmx_file("Cannot read energy file header. Corrupt file?");
 788         }
 789
 790         /* Now check whether this file is in single precision */
 791         if (!bWrongPrecision &&
 792             ((fr->e_size && (fr->nre == nre) &&
 793               (nre*4*(long int)sizeof(float) == fr->e_size)) ) )
 794         {
 795             fprintf(stderr, "Opened %s as single precision energy file\n", fn);
 796             free_enxnms(nre, nms);
 797         }
 798         else
 799         {
 800             gmx_fio_rewind(ef->fio);
 801             gmx_fio_checktype(ef->fio);
 802             gmx_fio_setprecision(ef->fio, TRUE);
 803             do_enxnms(ef, &nre, &nms);
 804             do_eheader(ef, &file_version, fr, nre, &bWrongPrecision, &bOK);
 805             if (!bOK)
 806             {
 807                 gmx_file("Cannot write energy file header; maybe you are out of disk space?");
 808             }
 809
 810             if (((fr->e_size && (fr->nre == nre) &&
 811                   (nre*4*(long int)sizeof(double) == fr->e_size)) ))
 812             {
 813                 fprintf(stderr, "Opened %s as double precision energy file\n",
 814                         fn);
 815             }
 816             else
 817             {
 818                 if (empty_file(fn))
 819                 {
 820                     gmx_fatal(FARGS, "File %s is empty", fn);
 821                 }
 822                 else
 823                 {
 824                     gmx_fatal(FARGS, "Energy file %s not recognized, maybe different CPU?",
 825                               fn);
 826                 }
 827             }
 828             free_enxnms(nre, nms);
 829         }
 830         free_enxframe(fr);
 831         sfree(fr);
 832         gmx_fio_rewind(ef->fio);
 833     }
 834     else
 835     {
 836         ef->fio = gmx_fio_open(fn, mode);
 837     }
 838
 839     ef->framenr   = 0;
 840     ef->frametime = 0;
 841     return ef;
 842 }
 843
 844 t_fileio *enx_file_pointer(const ener_file_t ef)
 845 {
 846     return ef->fio;
 847 }
 848
 849 static void convert_full_sums(ener_old_t *ener_old, t_enxframe *fr)
 850 {
 851     int    nstep_all;
 852     int    ne, ns, i;
 853     double esum_all, eav_all;
 854
 855     if (fr->nsum > 0)
 856     {
 857         ne = 0;
 858         ns = 0;
 859         for (i = 0; i < fr->nre; i++)
 860         {
 861             if (fr->ener[i].e    != 0)
 862             {
 863                 ne++;
 864             }
 865             if (fr->ener[i].esum != 0)
 866             {
 867                 ns++;
 868             }
 869         }
 870         if (ne > 0 && ns == 0)
 871         {
 872             /* We do not have all energy sums */
 873             fr->nsum = 0;
 874         }
 875     }
 876
 877     /* Convert old full simulation sums to sums between energy frames */
 878     nstep_all = fr->step - ener_old->first_step + 1;
 879     if (fr->nsum > 1 && fr->nsum == nstep_all && ener_old->nsum_prev > 0)
 880     {
 881         /* Set the new sum length: the frame step difference */
 882         fr->nsum = fr->step - ener_old->step_prev;
 883         for (i = 0; i < fr->nre; i++)
 884         {
 885             esum_all         = fr->ener[i].esum;
 886             eav_all          = fr->ener[i].eav;
 887             fr->ener[i].esum = esum_all - ener_old->ener_prev[i].esum;
 888             fr->ener[i].eav  = eav_all  - ener_old->ener_prev[i].eav
 889                 - dsqr(ener_old->ener_prev[i].esum/(nstep_all - fr->nsum)
 890                        - esum_all/nstep_all)*
 891                 (nstep_all - fr->nsum)*nstep_all/(double)fr->nsum;
 892             ener_old->ener_prev[i].esum = esum_all;
 893             ener_old->ener_prev[i].eav  = eav_all;
 894         }
 895         ener_old->nsum_prev = nstep_all;
 896     }
 897     else if (fr->nsum > 0)
 898     {
 899         if (fr->nsum != nstep_all)
 900         {
 901             fprintf(stderr, "\nWARNING: something is wrong with the energy sums, will not use exact averages\n");
 902             ener_old->nsum_prev = 0;
 903         }
 904         else
 905         {
 906             ener_old->nsum_prev = nstep_all;
 907         }
 908         /* Copy all sums to ener_prev */
 909         for (i = 0; i < fr->nre; i++)
 910         {
 911             ener_old->ener_prev[i].esum = fr->ener[i].esum;
 912             ener_old->ener_prev[i].eav  = fr->ener[i].eav;
 913         }
 914     }
 915
 916     ener_old->step_prev = fr->step;
 917 }
 918
 919 gmx_bool do_enx(ener_file_t ef, t_enxframe *fr)
 920 {
 921     int           file_version = -1;
 922     int           i, b;
 923     gmx_bool      bRead, bOK, bOK1, bSane;
 924     real          tmp1, tmp2, rdum;
 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=%"GMX_PRId64 ", nre=%d, nblock=%d, time=%g.\n"
 978                 "Trying to skip frame expect a crash though\n",
 979                 fr->step, 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 }