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