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/utility/futil.h"
  45 #include "gromacs/legacyheaders/macros.h"
  46
  47 #include "gromacs/fileio/gmxfio.h"
  48 #include "gromacs/fileio/xdrf.h"
  49 #include "gromacs/math/vec.h"
  50 #include "gromacs/pbcutil/pbc.h"
  51 #include "gromacs/topology/topology.h"
  52 #include "gromacs/utility/fatalerror.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     gmx_fio_checktype(ef->fio);
 392
 393     xdr = gmx_fio_getxdr(ef->fio);
 394
 395     if (!xdr_int(xdr, &magic))
 396     {
 397         if (!bRead)
 398         {
 399             gmx_file("Cannot write energy names to file; maybe you are out of disk space?");
 400         }
 401         *nre = 0;
 402         return;
 403     }
 404     if (magic > 0)
 405     {
 406         /* Assume this is an old edr format */
 407         file_version          = 1;
 408         *nre                  = magic;
 409         ef->eo.bOldFileOpen   = TRUE;
 410         ef->eo.bReadFirstStep = FALSE;
 411         srenew(ef->eo.ener_prev, *nre);
 412     }
 413     else
 414     {
 415         ef->eo.bOldFileOpen = FALSE;
 416
 417         if (magic != -55555)
 418         {
 419             gmx_fatal(FARGS, "Energy names magic number mismatch, this is not a GROMACS edr file");
 420         }
 421         file_version = enx_version;
 422         xdr_int(xdr, &file_version);
 423         if (file_version > enx_version)
 424         {
 425             gmx_fatal(FARGS, "reading tpx file (%s) version %d with version %d program", gmx_fio_getname(ef->fio), file_version, enx_version);
 426         }
 427         xdr_int(xdr, nre);
 428     }
 429     if (file_version != enx_version)
 430     {
 431         fprintf(stderr, "Note: enx file_version %d, software version %d\n",
 432                 file_version, enx_version);
 433     }
 434
 435     edr_strings(xdr, bRead, file_version, *nre, nms);
 436 }
 437
 438 static gmx_bool do_eheader(ener_file_t ef, int *file_version, t_enxframe *fr,
 439                            int nre_test, gmx_bool *bWrongPrecision, gmx_bool *bOK)
 440 {
 441     int          magic = -7777777;
 442     real         first_real_to_check;
 443     int          b, i, zero = 0, dum = 0;
 444     gmx_bool     bRead      = gmx_fio_getread(ef->fio);
 445     int          tempfix_nr = 0;
 446     int          ndisre     = 0;
 447     int          startb     = 0;
 448 #ifndef GMX_DOUBLE
 449     xdr_datatype dtreal = xdr_datatype_float;
 450 #else
 451     xdr_datatype dtreal = xdr_datatype_double;
 452 #endif
 453
 454     if (bWrongPrecision)
 455     {
 456         *bWrongPrecision = FALSE;
 457     }
 458
 459     *bOK = TRUE;
 460     /* The original energy frame started with a real,
 461      * so we have to use a real for compatibility.
 462      * This is VERY DIRTY code, since do_eheader can be called
 463      * with the wrong precision set and then we could read r > -1e10,
 464      * while actually the intention was r < -1e10.
 465      * When nre_test >= 0, do_eheader should therefore terminate
 466      * before the number of i/o calls starts depending on what has been read
 467      * (which is the case for for instance the block sizes for variable
 468      * number of blocks, where this number is read before).
 469      */
 470     first_real_to_check = -2e10;
 471     if (!gmx_fio_do_real(ef->fio, first_real_to_check))
 472     {
 473         return FALSE;
 474     }
 475     if (first_real_to_check > -1e10)
 476     {
 477         /* Assume we are reading an old format */
 478         *file_version = 1;
 479         fr->t         = first_real_to_check;
 480         if (!gmx_fio_do_int(ef->fio, dum))
 481         {
 482             *bOK = FALSE;
 483         }
 484         fr->step = dum;
 485     }
 486     else
 487     {
 488         if (!gmx_fio_do_int(ef->fio, magic))
 489         {
 490             *bOK = FALSE;
 491         }
 492         if (magic != -7777777)
 493         {
 494             enx_warning("Energy header magic number mismatch, this is not a GROMACS edr file");
 495             *bOK = FALSE;
 496             return FALSE;
 497         }
 498         *file_version = enx_version;
 499         if (!gmx_fio_do_int(ef->fio, *file_version))
 500         {
 501             *bOK = FALSE;
 502         }
 503         if (*bOK && *file_version > enx_version)
 504         {
 505             gmx_fatal(FARGS, "reading tpx file (%s) version %d with version %d program", gmx_fio_getname(ef->fio), file_version, enx_version);
 506         }
 507         if (!gmx_fio_do_double(ef->fio, fr->t))
 508         {
 509             *bOK = FALSE;
 510         }
 511         if (!gmx_fio_do_int64(ef->fio, fr->step))
 512         {
 513             *bOK = FALSE;
 514         }
 515         if (!bRead && fr->nsum == 1)
 516         {
 517             /* Do not store sums of length 1,
 518              * since this does not add information.
 519              */
 520             if (!gmx_fio_do_int(ef->fio, zero))
 521             {
 522                 *bOK = FALSE;
 523             }
 524         }
 525         else
 526         {
 527             if (!gmx_fio_do_int(ef->fio, fr->nsum))
 528             {
 529                 *bOK = FALSE;
 530             }
 531         }
 532         if (*file_version >= 3)
 533         {
 534             if (!gmx_fio_do_int64(ef->fio, fr->nsteps))
 535             {
 536                 *bOK = FALSE;
 537             }
 538         }
 539         else
 540         {
 541             fr->nsteps = max(1, fr->nsum);
 542         }
 543         if (*file_version >= 5)
 544         {
 545             if (!gmx_fio_do_double(ef->fio, fr->dt))
 546             {
 547                 *bOK = FALSE;
 548             }
 549         }
 550         else
 551         {
 552             fr->dt = 0;
 553         }
 554     }
 555     if (!gmx_fio_do_int(ef->fio, fr->nre))
 556     {
 557         *bOK = FALSE;
 558     }
 559     if (*file_version < 4)
 560     {
 561         if (!gmx_fio_do_int(ef->fio, ndisre))
 562         {
 563             *bOK = FALSE;
 564         }
 565     }
 566     else
 567     {
 568         /* now reserved for possible future use */
 569         if (!gmx_fio_do_int(ef->fio, dum))
 570         {
 571             *bOK = FALSE;
 572         }
 573     }
 574
 575     if (!gmx_fio_do_int(ef->fio, fr->nblock))
 576     {
 577         *bOK = FALSE;
 578     }
 579     if (fr->nblock < 0)
 580     {
 581         *bOK = FALSE;
 582     }
 583
 584     if (ndisre != 0)
 585     {
 586         if (*file_version >= 4)
 587         {
 588             enx_warning("Distance restraint blocks in old style in new style file");
 589             *bOK = FALSE;
 590             return FALSE;
 591         }
 592         fr->nblock += 1;
 593     }
 594
 595
 596     /* Frames could have nre=0, so we can not rely only on the fr->nre check */
 597     if (bRead && nre_test >= 0 &&
 598         ((fr->nre > 0 && fr->nre != nre_test) ||
 599          fr->nre < 0 || ndisre < 0 || fr->nblock < 0))
 600     {
 601         *bWrongPrecision = TRUE;
 602         return *bOK;
 603     }
 604
 605     /* we now know what these should be, or we've already bailed out because
 606        of wrong precision */
 607     if (*file_version == 1 && (fr->t < 0 || fr->t > 1e20 || fr->step < 0 ) )
 608     {
 609         enx_warning("edr file with negative step number or unreasonable time (and without version number).");
 610         *bOK = FALSE;
 611         return FALSE;
 612     }
 613
 614
 615     if (*bOK && bRead)
 616     {
 617         add_blocks_enxframe(fr, fr->nblock);
 618     }
 619
 620     startb = 0;
 621     if (ndisre > 0)
 622     {
 623         /* sub[0] is the instantaneous data, sub[1] is time averaged */
 624         add_subblocks_enxblock(&(fr->block[0]), 2);
 625         fr->block[0].id          = enxDISRE;
 626         fr->block[0].sub[0].nr   = ndisre;
 627         fr->block[0].sub[1].nr   = ndisre;
 628         fr->block[0].sub[0].type = dtreal;
 629         fr->block[0].sub[1].type = dtreal;
 630         startb++;
 631     }
 632
 633     /* read block header info */
 634     for (b = startb; b < fr->nblock; b++)
 635     {
 636         if (*file_version < 4)
 637         {
 638             /* blocks in old version files always have 1 subblock that
 639                consists of reals. */
 640             int nrint;
 641
 642             if (bRead)
 643             {
 644                 add_subblocks_enxblock(&(fr->block[b]), 1);
 645             }
 646             else
 647             {
 648                 if (fr->block[b].nsub != 1)
 649                 {
 650                     gmx_incons("Writing an old version .edr file with too many subblocks");
 651                 }
 652                 if (fr->block[b].sub[0].type != dtreal)
 653                 {
 654                     gmx_incons("Writing an old version .edr file the wrong subblock type");
 655                 }
 656             }
 657             nrint = fr->block[b].sub[0].nr;
 658
 659             if (!gmx_fio_do_int(ef->fio, nrint))
 660             {
 661                 *bOK = FALSE;
 662             }
 663             fr->block[b].id          = b - startb;
 664             fr->block[b].sub[0].nr   = nrint;
 665             fr->block[b].sub[0].type = dtreal;
 666         }
 667         else
 668         {
 669             int i;
 670             /* in the new version files, the block header only contains
 671                the ID and the number of subblocks */
 672             int nsub = fr->block[b].nsub;
 673             *bOK = *bOK && gmx_fio_do_int(ef->fio, fr->block[b].id);
 674             *bOK = *bOK && gmx_fio_do_int(ef->fio, nsub);
 675
 676             fr->block[b].nsub = nsub;
 677             if (bRead)
 678             {
 679                 add_subblocks_enxblock(&(fr->block[b]), nsub);
 680             }
 681
 682             /* read/write type & size for each subblock */
 683             for (i = 0; i < nsub; i++)
 684             {
 685                 t_enxsubblock *sub    = &(fr->block[b].sub[i]); /* shortcut */
 686                 int            typenr = sub->type;
 687
 688                 *bOK = *bOK && gmx_fio_do_int(ef->fio, typenr);
 689                 *bOK = *bOK && gmx_fio_do_int(ef->fio, sub->nr);
 690
 691                 sub->type = (xdr_datatype)typenr;
 692             }
 693         }
 694     }
 695     if (!gmx_fio_do_int(ef->fio, fr->e_size))
 696     {
 697         *bOK = FALSE;
 698     }
 699
 700     /* now reserved for possible future use */
 701     if (!gmx_fio_do_int(ef->fio, dum))
 702     {
 703         *bOK = FALSE;
 704     }
 705
 706     /* Do a dummy int to keep the format compatible with the old code */
 707     if (!gmx_fio_do_int(ef->fio, dum))
 708     {
 709         *bOK = FALSE;
 710     }
 711
 712     if (*bOK && *file_version == 1 && nre_test < 0)
 713     {
 714         if (!ef->eo.bReadFirstStep)
 715         {
 716             ef->eo.bReadFirstStep = TRUE;
 717             ef->eo.first_step     = fr->step;
 718             ef->eo.step_prev      = fr->step;
 719             ef->eo.nsum_prev      = 0;
 720         }
 721
 722         fr->nsum   = fr->step - ef->eo.first_step + 1;
 723         fr->nsteps = fr->step - ef->eo.step_prev;
 724         fr->dt     = 0;
 725     }
 726
 727     return *bOK;
 728 }
 729
 730 void free_enxnms(int n, gmx_enxnm_t *nms)
 731 {
 732     int i;
 733
 734     for (i = 0; i < n; i++)
 735     {
 736         sfree(nms[i].name);
 737         sfree(nms[i].unit);
 738     }
 739
 740     sfree(nms);
 741 }
 742
 743 void close_enx(ener_file_t ef)
 744 {
 745     if (gmx_fio_close(ef->fio) != 0)
 746     {
 747         gmx_file("Cannot close energy file; it might be corrupt, or maybe you are out of disk space?");
 748     }
 749 }
 750
 751 static gmx_bool empty_file(const char *fn)
 752 {
 753     FILE    *fp;
 754     char     dum;
 755     int      ret;
 756     gmx_bool bEmpty;
 757
 758     fp     = gmx_fio_fopen(fn, "r");
 759     ret    = fread(&dum, sizeof(dum), 1, fp);
 760     bEmpty = feof(fp);
 761     gmx_fio_fclose(fp);
 762
 763     return bEmpty;
 764 }
 765
 766
 767 ener_file_t open_enx(const char *fn, const char *mode)
 768 {
 769     int               nre, i;
 770     gmx_enxnm_t      *nms          = NULL;
 771     int               file_version = -1;
 772     t_enxframe       *fr;
 773     gmx_bool          bWrongPrecision, bOK = TRUE;
 774     struct ener_file *ef;
 775
 776     snew(ef, 1);
 777
 778     if (mode[0] == 'r')
 779     {
 780         ef->fio = gmx_fio_open(fn, mode);
 781         gmx_fio_checktype(ef->fio);
 782         gmx_fio_setprecision(ef->fio, FALSE);
 783         do_enxnms(ef, &nre, &nms);
 784         snew(fr, 1);
 785         do_eheader(ef, &file_version, fr, nre, &bWrongPrecision, &bOK);
 786         if (!bOK)
 787         {
 788             gmx_file("Cannot read energy file header. Corrupt file?");
 789         }
 790
 791         /* Now check whether this file is in single precision */
 792         if (!bWrongPrecision &&
 793             ((fr->e_size && (fr->nre == nre) &&
 794               (nre*4*(long int)sizeof(float) == fr->e_size)) ) )
 795         {
 796             fprintf(stderr, "Opened %s as single precision energy file\n", fn);
 797             free_enxnms(nre, nms);
 798         }
 799         else
 800         {
 801             gmx_fio_rewind(ef->fio);
 802             gmx_fio_checktype(ef->fio);
 803             gmx_fio_setprecision(ef->fio, TRUE);
 804             do_enxnms(ef, &nre, &nms);
 805             do_eheader(ef, &file_version, fr, nre, &bWrongPrecision, &bOK);
 806             if (!bOK)
 807             {
 808                 gmx_file("Cannot write energy file header; maybe you are out of disk space?");
 809             }
 810
 811             if (((fr->e_size && (fr->nre == nre) &&
 812                   (nre*4*(long int)sizeof(double) == fr->e_size)) ))
 813             {
 814                 fprintf(stderr, "Opened %s as double precision energy file\n",
 815                         fn);
 816             }
 817             else
 818             {
 819                 if (empty_file(fn))
 820                 {
 821                     gmx_fatal(FARGS, "File %s is empty", fn);
 822                 }
 823                 else
 824                 {
 825                     gmx_fatal(FARGS, "Energy file %s not recognized, maybe different CPU?",
 826                               fn);
 827                 }
 828             }
 829             free_enxnms(nre, nms);
 830         }
 831         free_enxframe(fr);
 832         sfree(fr);
 833         gmx_fio_rewind(ef->fio);
 834     }
 835     else
 836     {
 837         ef->fio = gmx_fio_open(fn, mode);
 838     }
 839
 840     ef->framenr   = 0;
 841     ef->frametime = 0;
 842     return ef;
 843 }
 844
 845 t_fileio *enx_file_pointer(const ener_file_t ef)
 846 {
 847     return ef->fio;
 848 }
 849
 850 static void convert_full_sums(ener_old_t *ener_old, t_enxframe *fr)
 851 {
 852     int    nstep_all;
 853     int    ne, ns, i;
 854     double esum_all, eav_all;
 855
 856     if (fr->nsum > 0)
 857     {
 858         ne = 0;
 859         ns = 0;
 860         for (i = 0; i < fr->nre; i++)
 861         {
 862             if (fr->ener[i].e    != 0)
 863             {
 864                 ne++;
 865             }
 866             if (fr->ener[i].esum != 0)
 867             {
 868                 ns++;
 869             }
 870         }
 871         if (ne > 0 && ns == 0)
 872         {
 873             /* We do not have all energy sums */
 874             fr->nsum = 0;
 875         }
 876     }
 877
 878     /* Convert old full simulation sums to sums between energy frames */
 879     nstep_all = fr->step - ener_old->first_step + 1;
 880     if (fr->nsum > 1 && fr->nsum == nstep_all && ener_old->nsum_prev > 0)
 881     {
 882         /* Set the new sum length: the frame step difference */
 883         fr->nsum = fr->step - ener_old->step_prev;
 884         for (i = 0; i < fr->nre; i++)
 885         {
 886             esum_all         = fr->ener[i].esum;
 887             eav_all          = fr->ener[i].eav;
 888             fr->ener[i].esum = esum_all - ener_old->ener_prev[i].esum;
 889             fr->ener[i].eav  = eav_all  - ener_old->ener_prev[i].eav
 890                 - dsqr(ener_old->ener_prev[i].esum/(nstep_all - fr->nsum)
 891                        - esum_all/nstep_all)*
 892                 (nstep_all - fr->nsum)*nstep_all/(double)fr->nsum;
 893             ener_old->ener_prev[i].esum = esum_all;
 894             ener_old->ener_prev[i].eav  = eav_all;
 895         }
 896         ener_old->nsum_prev = nstep_all;
 897     }
 898     else if (fr->nsum > 0)
 899     {
 900         if (fr->nsum != nstep_all)
 901         {
 902             fprintf(stderr, "\nWARNING: something is wrong with the energy sums, will not use exact averages\n");
 903             ener_old->nsum_prev = 0;
 904         }
 905         else
 906         {
 907             ener_old->nsum_prev = nstep_all;
 908         }
 909         /* Copy all sums to ener_prev */
 910         for (i = 0; i < fr->nre; i++)
 911         {
 912             ener_old->ener_prev[i].esum = fr->ener[i].esum;
 913             ener_old->ener_prev[i].eav  = fr->ener[i].eav;
 914         }
 915     }
 916
 917     ener_old->step_prev = fr->step;
 918 }
 919
 920 gmx_bool do_enx(ener_file_t ef, t_enxframe *fr)
 921 {
 922     int           file_version = -1;
 923     int           i, b;
 924     gmx_bool      bRead, bOK, bOK1, bSane;
 925     real          tmp1, tmp2, rdum;
 926     /*int       d_size;*/
 927
 928     bOK   = TRUE;
 929     bRead = gmx_fio_getread(ef->fio);
 930     if (!bRead)
 931     {
 932         fr->e_size = fr->nre*sizeof(fr->ener[0].e)*4;
 933         /*d_size = fr->ndisre*(sizeof(real)*2);*/
 934     }
 935     gmx_fio_checktype(ef->fio);
 936
 937     if (!do_eheader(ef, &file_version, fr, -1, NULL, &bOK))
 938     {
 939         if (bRead)
 940         {
 941             fprintf(stderr, "\rLast energy frame read %d time %8.3f         ",
 942                     ef->framenr-1, ef->frametime);
 943             if (!bOK)
 944             {
 945                 fprintf(stderr,
 946                         "\nWARNING: Incomplete energy frame: nr %d time %8.3f\n",
 947                         ef->framenr, fr->t);
 948             }
 949         }
 950         else
 951         {
 952             gmx_file("Cannot write energy file header; maybe you are out of disk space?");
 953         }
 954         return FALSE;
 955     }
 956     if (bRead)
 957     {
 958         if ((ef->framenr <   20 || ef->framenr %   10 == 0) &&
 959             (ef->framenr <  200 || ef->framenr %  100 == 0) &&
 960             (ef->framenr < 2000 || ef->framenr % 1000 == 0))
 961         {
 962             fprintf(stderr, "\rReading energy frame %6d time %8.3f         ",
 963                     ef->framenr, fr->t);
 964         }
 965         ef->framenr++;
 966         ef->frametime = fr->t;
 967     }
 968     /* Check sanity of this header */
 969     bSane = fr->nre > 0;
 970     for (b = 0; b < fr->nblock; b++)
 971     {
 972         bSane = bSane || (fr->block[b].nsub > 0);
 973     }
 974     if (!((fr->step >= 0) && bSane))
 975     {
 976         fprintf(stderr, "\nWARNING: there may be something wrong with energy file %s\n",
 977                 gmx_fio_getname(ef->fio));
 978         fprintf(stderr, "Found: step=%"GMX_PRId64 ", nre=%d, nblock=%d, time=%g.\n"
 979                 "Trying to skip frame expect a crash though\n",
 980                 fr->step, fr->nre, fr->nblock, fr->t);
 981     }
 982     if (bRead && fr->nre > fr->e_alloc)
 983     {
 984         srenew(fr->ener, fr->nre);
 985         for (i = fr->e_alloc; (i < fr->nre); i++)
 986         {
 987             fr->ener[i].e    = 0;
 988             fr->ener[i].eav  = 0;
 989             fr->ener[i].esum = 0;
 990         }
 991         fr->e_alloc = fr->nre;
 992     }
 993
 994     for (i = 0; i < fr->nre; i++)
 995     {
 996         bOK = bOK && gmx_fio_do_real(ef->fio, fr->ener[i].e);
 997
 998         /* Do not store sums of length 1,
 999          * since this does not add information.
1000          */
1001         if (file_version == 1 ||
1002             (bRead && fr->nsum > 0) || fr->nsum > 1)
1003         {
1004             tmp1 = fr->ener[i].eav;
1005             bOK  = bOK && gmx_fio_do_real(ef->fio, tmp1);
1006             if (bRead)
1007             {
1008                 fr->ener[i].eav = tmp1;
1009             }
1010
1011             /* This is to save only in single precision (unless compiled in DP) */
1012             tmp2 = fr->ener[i].esum;
1013             bOK  = bOK && gmx_fio_do_real(ef->fio, tmp2);
1014             if (bRead)
1015             {
1016                 fr->ener[i].esum = tmp2;
1017             }
1018
1019             if (file_version == 1)
1020             {
1021                 /* Old, unused real */
1022                 rdum = 0;
1023                 bOK  = bOK && gmx_fio_do_real(ef->fio, rdum);
1024             }
1025         }
1026     }
1027
1028     /* Here we can not check for file_version==1, since one could have
1029      * continued an old format simulation with a new one with mdrun -append.
1030      */
1031     if (bRead && ef->eo.bOldFileOpen)
1032     {
1033         /* Convert old full simulation sums to sums between energy frames */
1034         convert_full_sums(&(ef->eo), fr);
1035     }
1036     /* read the blocks */
1037     for (b = 0; b < fr->nblock; b++)
1038     {
1039         /* now read the subblocks. */
1040         int nsub = fr->block[b].nsub; /* shortcut */
1041         int i;
1042
1043         for (i = 0; i < nsub; i++)
1044         {
1045             t_enxsubblock *sub = &(fr->block[b].sub[i]); /* shortcut */
1046
1047             if (bRead)
1048             {
1049                 enxsubblock_alloc(sub);
1050             }
1051
1052             /* read/write data */
1053             bOK1 = TRUE;
1054             switch (sub->type)
1055             {
1056                 case xdr_datatype_float:
1057                     bOK1 = gmx_fio_ndo_float(ef->fio, sub->fval, sub->nr);
1058                     break;
1059                 case xdr_datatype_double:
1060                     bOK1 = gmx_fio_ndo_double(ef->fio, sub->dval, sub->nr);
1061                     break;
1062                 case xdr_datatype_int:
1063                     bOK1 = gmx_fio_ndo_int(ef->fio, sub->ival, sub->nr);
1064                     break;
1065                 case xdr_datatype_int64:
1066                     bOK1 = gmx_fio_ndo_int64(ef->fio, sub->lval, sub->nr);
1067                     break;
1068                 case xdr_datatype_char:
1069                     bOK1 = gmx_fio_ndo_uchar(ef->fio, sub->cval, sub->nr);
1070                     break;
1071                 case xdr_datatype_string:
1072                     bOK1 = gmx_fio_ndo_string(ef->fio, sub->sval, sub->nr);
1073                     break;
1074                 default:
1075                     gmx_incons("Reading unknown block data type: this file is corrupted or from the future");
1076             }
1077             bOK = bOK && bOK1;
1078         }
1079     }
1080
1081     if (!bRead)
1082     {
1083         if (gmx_fio_flush(ef->fio) != 0)
1084         {
1085             gmx_file("Cannot write energy file; maybe you are out of disk space?");
1086         }
1087     }
1088
1089     if (!bOK)
1090     {
1091         if (bRead)
1092         {
1093             fprintf(stderr, "\nLast energy frame read %d",
1094                     ef->framenr-1);
1095             fprintf(stderr, "\nWARNING: Incomplete energy frame: nr %d time %8.3f\n",
1096                     ef->framenr, fr->t);
1097         }
1098         else
1099         {
1100             gmx_fatal(FARGS, "could not write energies");
1101         }
1102         return FALSE;
1103     }
1104
1105     return TRUE;
1106 }
1107
1108 static real find_energy(const char *name, int nre, gmx_enxnm_t *enm,
1109                         t_enxframe *fr)
1110 {
1111     int i;
1112
1113     for (i = 0; i < nre; i++)
1114     {
1115         if (strcmp(enm[i].name, name) == 0)
1116         {
1117             return fr->ener[i].e;
1118         }
1119     }
1120
1121     gmx_fatal(FARGS, "Could not find energy term named '%s'", name);
1122
1123     return 0;
1124 }
1125
1126
1127 void get_enx_state(const char *fn, real t, gmx_groups_t *groups, t_inputrec *ir,
1128                    t_state *state)
1129 {
1130     /* Should match the names in mdebin.c */
1131     static const char *boxvel_nm[] = {
1132         "Box-Vel-XX", "Box-Vel-YY", "Box-Vel-ZZ",
1133         "Box-Vel-YX", "Box-Vel-ZX", "Box-Vel-ZY"
1134     };
1135
1136     static const char *pcouplmu_nm[] = {
1137         "Pcoupl-Mu-XX", "Pcoupl-Mu-YY", "Pcoupl-Mu-ZZ",
1138         "Pcoupl-Mu-YX", "Pcoupl-Mu-ZX", "Pcoupl-Mu-ZY"
1139     };
1140     static const char *baro_nm[] = {
1141         "Barostat"
1142     };
1143
1144
1145     int          ind0[] = { XX, YY, ZZ, YY, ZZ, ZZ };
1146     int          ind1[] = { XX, YY, ZZ, XX, XX, YY };
1147     int          nre, nfr, i, j, ni, npcoupl;
1148     char         buf[STRLEN];
1149     const char  *bufi;
1150     gmx_enxnm_t *enm = NULL;
1151     t_enxframe  *fr;
1152     ener_file_t  in;
1153
1154     in = open_enx(fn, "r");
1155     do_enxnms(in, &nre, &enm);
1156     snew(fr, 1);
1157     nfr = 0;
1158     while ((nfr == 0 || fr->t != t) && do_enx(in, fr))
1159     {
1160         nfr++;
1161     }
1162     close_enx(in);
1163     fprintf(stderr, "\n");
1164
1165     if (nfr == 0 || fr->t != t)
1166     {
1167         gmx_fatal(FARGS, "Could not find frame with time %f in '%s'", t, fn);
1168     }
1169
1170     npcoupl = TRICLINIC(ir->compress) ? 6 : 3;
1171     if (ir->epc == epcPARRINELLORAHMAN)
1172     {
1173         clear_mat(state->boxv);
1174         for (i = 0; i < npcoupl; i++)
1175         {
1176             state->boxv[ind0[i]][ind1[i]] =
1177                 find_energy(boxvel_nm[i], nre, enm, fr);
1178         }
1179         fprintf(stderr, "\nREAD %d BOX VELOCITIES FROM %s\n\n", npcoupl, fn);
1180     }
1181
1182     if (ir->etc == etcNOSEHOOVER)
1183     {
1184         char cns[20];
1185
1186         cns[0] = '\0';
1187
1188         for (i = 0; i < state->ngtc; i++)
1189         {
1190             ni   = groups->grps[egcTC].nm_ind[i];
1191             bufi = *(groups->grpname[ni]);
1192             for (j = 0; (j < state->nhchainlength); j++)
1193             {
1194                 if (IR_NVT_TROTTER(ir))
1195                 {
1196                     sprintf(cns, "-%d", j);
1197                 }
1198                 sprintf(buf, "Xi%s-%s", cns, bufi);
1199                 state->nosehoover_xi[i] = find_energy(buf, nre, enm, fr);
1200                 sprintf(buf, "vXi%s-%s", cns, bufi);
1201                 state->nosehoover_vxi[i] = find_energy(buf, nre, enm, fr);
1202             }
1203
1204         }
1205         fprintf(stderr, "\nREAD %d NOSE-HOOVER Xi chains FROM %s\n\n", state->ngtc, fn);
1206
1207         if (IR_NPT_TROTTER(ir) || IR_NPH_TROTTER(ir))
1208         {
1209             for (i = 0; i < state->nnhpres; i++)
1210             {
1211                 bufi = baro_nm[0]; /* All barostat DOF's together for now */
1212                 for (j = 0; (j < state->nhchainlength); j++)
1213                 {
1214                     sprintf(buf, "Xi-%d-%s", j, bufi);
1215                     state->nhpres_xi[i] = find_energy(buf, nre, enm, fr);
1216                     sprintf(buf, "vXi-%d-%s", j, bufi);
1217                     state->nhpres_vxi[i] = find_energy(buf, nre, enm, fr);
1218                 }
1219             }
1220             fprintf(stderr, "\nREAD %d NOSE-HOOVER BAROSTAT Xi chains FROM %s\n\n", state->nnhpres, fn);
1221         }
1222     }
1223
1224     free_enxnms(nre, enm);
1225     free_enxframe(fr);
1226     sfree(fr);
1227 }