src/gromacs/fileio/enxio.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include "enxio.h"
  40
  41 #include <cstdlib>
  42 #include <cstring>
  43
  44 #include <algorithm>
  45
  46 #include "gromacs/fileio/gmxfio.h"
  47 #include "gromacs/fileio/gmxfio-xdr.h"
  48 #include "gromacs/fileio/xdrf.h"
  49 #include "gromacs/legacyheaders/macros.h"
  50 #include "gromacs/math/vec.h"
  51 #include "gromacs/pbcutil/pbc.h"
  52 #include "gromacs/topology/topology.h"
  53 #include "gromacs/utility/fatalerror.h"
  54 #include "gromacs/utility/futil.h"
  55 #include "gromacs/utility/gmxassert.h"
  56 #include "gromacs/utility/smalloc.h"
  57
  58 /* The source code in this file should be thread-safe.
  59          Please keep it that way. */
  60
  61 /* This number should be increased whenever the file format changes! */
  62 static const int enx_version = 5;
  63
  64 const char      *enx_block_id_name[] = {
  65     "Averaged orientation restraints",
  66     "Instantaneous orientation restraints",
  67     "Orientation restraint order tensor(s)",
  68     "Distance restraints",
  69     "Free energy data",
  70     "BAR histogram",
  71     "Delta H raw data"
  72 };
  73
  74
  75 /* Stuff for reading pre 4.1 energy files */
  76 typedef struct {
  77     gmx_bool     bOldFileOpen;   /* Is this an open old file? */
  78     gmx_bool     bReadFirstStep; /* Did we read the first step? */
  79     int          first_step;     /* First step in the energy file */
  80     int          step_prev;      /* Previous step */
  81     int          nsum_prev;      /* Previous step sum length */
  82     t_energy    *ener_prev;      /* Previous energy sums */
  83 } ener_old_t;
  84
  85 struct ener_file
  86 {
  87     ener_old_t eo;
  88     t_fileio  *fio;
  89     int        framenr;
  90     real       frametime;
  91 };
  92
  93 static void enxsubblock_init(t_enxsubblock *sb)
  94 {
  95     sb->nr = 0;
  96 #ifdef GMX_DOUBLE
  97     sb->type = xdr_datatype_double;
  98 #else
  99     sb->type = xdr_datatype_float;
 100 #endif
 101     sb->fval       = NULL;
 102     sb->dval       = NULL;
 103     sb->ival       = NULL;
 104     sb->lval       = NULL;
 105     sb->cval       = NULL;
 106     sb->sval       = NULL;
 107     sb->fval_alloc = 0;
 108     sb->dval_alloc = 0;
 109     sb->ival_alloc = 0;
 110     sb->lval_alloc = 0;
 111     sb->cval_alloc = 0;
 112     sb->sval_alloc = 0;
 113 }
 114
 115 static void enxsubblock_free(t_enxsubblock *sb)
 116 {
 117     if (sb->fval_alloc)
 118     {
 119         sfree(sb->fval);
 120         sb->fval_alloc = 0;
 121         sb->fval       = NULL;
 122     }
 123     if (sb->dval_alloc)
 124     {
 125         sfree(sb->dval);
 126         sb->dval_alloc = 0;
 127         sb->dval       = NULL;
 128     }
 129     if (sb->ival_alloc)
 130     {
 131         sfree(sb->ival);
 132         sb->ival_alloc = 0;
 133         sb->ival       = NULL;
 134     }
 135     if (sb->lval_alloc)
 136     {
 137         sfree(sb->lval);
 138         sb->lval_alloc = 0;
 139         sb->lval       = NULL;
 140     }
 141     if (sb->cval_alloc)
 142     {
 143         sfree(sb->cval);
 144         sb->cval_alloc = 0;
 145         sb->cval       = NULL;
 146     }
 147     if (sb->sval_alloc)
 148     {
 149         int i;
 150
 151         for (i = 0; i < sb->sval_alloc; i++)
 152         {
 153             if (sb->sval[i])
 154             {
 155                 sfree(sb->sval[i]);
 156             }
 157         }
 158         sfree(sb->sval);
 159         sb->sval_alloc = 0;
 160         sb->sval       = NULL;
 161     }
 162 }
 163
 164 /* allocate the appropriate amount of memory for the given type and nr */
 165 static void enxsubblock_alloc(t_enxsubblock *sb)
 166 {
 167     /* allocate the appropriate amount of memory */
 168     switch (sb->type)
 169     {
 170         case xdr_datatype_float:
 171             if (sb->nr > sb->fval_alloc)
 172             {
 173                 srenew(sb->fval, sb->nr);
 174                 sb->fval_alloc = sb->nr;
 175             }
 176             break;
 177         case xdr_datatype_double:
 178             if (sb->nr > sb->dval_alloc)
 179             {
 180                 srenew(sb->dval, sb->nr);
 181                 sb->dval_alloc = sb->nr;
 182             }
 183             break;
 184         case xdr_datatype_int:
 185             if (sb->nr > sb->ival_alloc)
 186             {
 187                 srenew(sb->ival, sb->nr);
 188                 sb->ival_alloc = sb->nr;
 189             }
 190             break;
 191         case xdr_datatype_int64:
 192             if (sb->nr > sb->lval_alloc)
 193             {
 194                 srenew(sb->lval, sb->nr);
 195                 sb->lval_alloc = sb->nr;
 196             }
 197             break;
 198         case xdr_datatype_char:
 199             if (sb->nr > sb->cval_alloc)
 200             {
 201                 srenew(sb->cval, sb->nr);
 202                 sb->cval_alloc = sb->nr;
 203             }
 204             break;
 205         case xdr_datatype_string:
 206             if (sb->nr > sb->sval_alloc)
 207             {
 208                 int i;
 209
 210                 srenew(sb->sval, sb->nr);
 211                 for (i = sb->sval_alloc; i < sb->nr; i++)
 212                 {
 213                     sb->sval[i] = NULL;
 214                 }
 215                 sb->sval_alloc = sb->nr;
 216             }
 217             break;
 218         default:
 219             gmx_incons("Unknown block type: this file is corrupted or from the future");
 220     }
 221 }
 222
 223 static void enxblock_init(t_enxblock *eb)
 224 {
 225     eb->id         = enxOR;
 226     eb->nsub       = 0;
 227     eb->sub        = NULL;
 228     eb->nsub_alloc = 0;
 229 }
 230
 231 static void enxblock_free(t_enxblock *eb)
 232 {
 233     if (eb->nsub_alloc > 0)
 234     {
 235         int i;
 236         for (i = 0; i < eb->nsub_alloc; i++)
 237         {
 238             enxsubblock_free(&(eb->sub[i]));
 239         }
 240         sfree(eb->sub);
 241         eb->nsub_alloc = 0;
 242         eb->sub        = NULL;
 243     }
 244 }
 245
 246 void init_enxframe(t_enxframe *fr)
 247 {
 248     fr->e_alloc = 0;
 249     fr->ener    = NULL;
 250
 251     /*fr->d_alloc=0;*/
 252     fr->ener = NULL;
 253
 254     /*fr->ndisre=0;*/
 255
 256     fr->nblock       = 0;
 257     fr->nblock_alloc = 0;
 258     fr->block        = NULL;
 259 }
 260
 261
 262 void free_enxframe(t_enxframe *fr)
 263 {
 264     int b;
 265
 266     if (fr->e_alloc)
 267     {
 268         sfree(fr->ener);
 269     }
 270     for (b = 0; b < fr->nblock_alloc; b++)
 271     {
 272         enxblock_free(&(fr->block[b]));
 273     }
 274     sfree(fr->block);
 275 }
 276
 277 void add_blocks_enxframe(t_enxframe *fr, int n)
 278 {
 279     fr->nblock = n;
 280     if (n > fr->nblock_alloc)
 281     {
 282         int b;
 283
 284         srenew(fr->block, n);
 285         for (b = fr->nblock_alloc; b < fr->nblock; b++)
 286         {
 287             enxblock_init(&(fr->block[b]));
 288         }
 289         fr->nblock_alloc = n;
 290     }
 291 }
 292
 293 t_enxblock *find_block_id_enxframe(t_enxframe *ef, int id, t_enxblock *prev)
 294 {
 295     gmx_off_t starti = 0;
 296     gmx_off_t i;
 297
 298     if (prev)
 299     {
 300         starti = (prev - ef->block) + 1;
 301     }
 302     for (i = starti; i < ef->nblock; i++)
 303     {
 304         if (ef->block[i].id == id)
 305         {
 306             return &(ef->block[i]);
 307         }
 308     }
 309     return NULL;
 310 }
 311
 312 void add_subblocks_enxblock(t_enxblock *eb, int n)
 313 {
 314     eb->nsub = n;
 315     if (eb->nsub > eb->nsub_alloc)
 316     {
 317         int b;
 318
 319         srenew(eb->sub, n);
 320         for (b = eb->nsub_alloc; b < n; b++)
 321         {
 322             enxsubblock_init(&(eb->sub[b]));
 323         }
 324         eb->nsub_alloc = n;
 325     }
 326 }
 327
 328 static void enx_warning(const char *msg)
 329 {
 330     if (getenv("GMX_ENX_NO_FATAL") != NULL)
 331     {
 332         gmx_warning(msg);
 333     }
 334     else
 335     {
 336         gmx_fatal(FARGS, "%s\n%s",
 337                   msg,
 338                   "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");
 339     }
 340 }
 341
 342 static void edr_strings(XDR *xdr, gmx_bool bRead, int file_version,
 343                         int n, gmx_enxnm_t **nms)
 344 {
 345     int          i;
 346     gmx_enxnm_t *nm;
 347
 348     if (*nms == NULL)
 349     {
 350         snew(*nms, n);
 351     }
 352     for (i = 0; i < n; i++)
 353     {
 354         nm = &(*nms)[i];
 355         if (bRead)
 356         {
 357             if (nm->name)
 358             {
 359                 sfree(nm->name);
 360                 nm->name = NULL;
 361             }
 362             if (nm->unit)
 363             {
 364                 sfree(nm->unit);
 365                 nm->unit = NULL;
 366             }
 367         }
 368         if (!xdr_string(xdr, &(nm->name), STRLEN))
 369         {
 370             gmx_file("Cannot write energy names to file; maybe you are out of disk space?");
 371         }
 372         if (file_version >= 2)
 373         {
 374             if (!xdr_string(xdr, &(nm->unit), STRLEN))
 375             {
 376                 gmx_file("Cannot write energy names to file; maybe you are out of disk space?");
 377             }
 378         }
 379         else
 380         {
 381             nm->unit = gmx_strdup("kJ/mol");
 382         }
 383     }
 384 }
 385
 386 void do_enxnms(ener_file_t ef, int *nre, gmx_enxnm_t **nms)
 387 {
 388     int      magic = -55555;
 389     XDR     *xdr;
 390     gmx_bool bRead = gmx_fio_getread(ef->fio);
 391     int      file_version;
 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, zero = 0, dum = 0;
 444     gmx_bool     bRead      = gmx_fio_getread(ef->fio);
 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 = std::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         if (bWrongPrecision)
 601         {
 602             *bWrongPrecision = TRUE;
 603         }
 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 /*!\brief Return TRUE if a file exists but is empty, otherwise FALSE.
 754  *
 755  * If the file exists but has length larger than zero, if it does not exist, or
 756  * if there is a file system error, FALSE will be returned instead.
 757  *
 758  * \param fn  File name to test
 759  *
 760  * \return TRUE if file could be open but is empty, otherwise FALSE.
 761  */
 762 static gmx_bool
 763 empty_file(const char *fn)
 764 {
 765     FILE    *fp;
 766     char     dum;
 767     int      ret;
 768     gmx_bool bEmpty;
 769
 770     fp     = gmx_fio_fopen(fn, "r");
 771     ret    = fread(&dum, sizeof(dum), 1, fp);
 772     bEmpty = feof(fp);
 773     gmx_fio_fclose(fp);
 774
 775     // bEmpty==TRUE but ret!=0 would likely be some strange I/O error, but at
 776     // least it's not a normal empty file, so we return FALSE in that case.
 777     return (bEmpty && ret == 0);
 778 }
 779
 780
 781 ener_file_t open_enx(const char *fn, const char *mode)
 782 {
 783     int               nre;
 784     gmx_enxnm_t      *nms          = NULL;
 785     int               file_version = -1;
 786     t_enxframe       *fr;
 787     gmx_bool          bWrongPrecision, bOK = TRUE;
 788     struct ener_file *ef;
 789
 790     snew(ef, 1);
 791
 792     if (mode[0] == 'r')
 793     {
 794         ef->fio = gmx_fio_open(fn, mode);
 795         gmx_fio_setprecision(ef->fio, FALSE);
 796         do_enxnms(ef, &nre, &nms);
 797         snew(fr, 1);
 798         do_eheader(ef, &file_version, fr, nre, &bWrongPrecision, &bOK);
 799         if (!bOK)
 800         {
 801             gmx_file("Cannot read energy file header. Corrupt file?");
 802         }
 803
 804         /* Now check whether this file is in single precision */
 805         if (!bWrongPrecision &&
 806             ((fr->e_size && (fr->nre == nre) &&
 807               (nre*4*(long int)sizeof(float) == fr->e_size)) ) )
 808         {
 809             fprintf(stderr, "Opened %s as single precision energy file\n", fn);
 810             free_enxnms(nre, nms);
 811         }
 812         else
 813         {
 814             gmx_fio_rewind(ef->fio);
 815             gmx_fio_setprecision(ef->fio, TRUE);
 816             do_enxnms(ef, &nre, &nms);
 817             do_eheader(ef, &file_version, fr, nre, &bWrongPrecision, &bOK);
 818             if (!bOK)
 819             {
 820                 gmx_file("Cannot write energy file header; maybe you are out of disk space?");
 821             }
 822
 823             if (((fr->e_size && (fr->nre == nre) &&
 824                   (nre*4*(long int)sizeof(double) == fr->e_size)) ))
 825             {
 826                 fprintf(stderr, "Opened %s as double precision energy file\n",
 827                         fn);
 828             }
 829             else
 830             {
 831                 if (empty_file(fn))
 832                 {
 833                     gmx_fatal(FARGS, "File %s is empty", fn);
 834                 }
 835                 else
 836                 {
 837                     gmx_fatal(FARGS, "Energy file %s not recognized, maybe different CPU?",
 838                               fn);
 839                 }
 840             }
 841             free_enxnms(nre, nms);
 842         }
 843         free_enxframe(fr);
 844         sfree(fr);
 845         gmx_fio_rewind(ef->fio);
 846     }
 847     else
 848     {
 849         ef->fio = gmx_fio_open(fn, mode);
 850     }
 851
 852     ef->framenr   = 0;
 853     ef->frametime = 0;
 854     return ef;
 855 }
 856
 857 t_fileio *enx_file_pointer(const ener_file_t ef)
 858 {
 859     return ef->fio;
 860 }
 861
 862 static void convert_full_sums(ener_old_t *ener_old, t_enxframe *fr)
 863 {
 864     int    nstep_all;
 865     int    ne, ns, i;
 866     double esum_all, eav_all;
 867
 868     if (fr->nsum > 0)
 869     {
 870         ne = 0;
 871         ns = 0;
 872         for (i = 0; i < fr->nre; i++)
 873         {
 874             if (fr->ener[i].e    != 0)
 875             {
 876                 ne++;
 877             }
 878             if (fr->ener[i].esum != 0)
 879             {
 880                 ns++;
 881             }
 882         }
 883         if (ne > 0 && ns == 0)
 884         {
 885             /* We do not have all energy sums */
 886             fr->nsum = 0;
 887         }
 888     }
 889
 890     /* Convert old full simulation sums to sums between energy frames */
 891     nstep_all = fr->step - ener_old->first_step + 1;
 892     if (fr->nsum > 1 && fr->nsum == nstep_all && ener_old->nsum_prev > 0)
 893     {
 894         /* Set the new sum length: the frame step difference */
 895         fr->nsum = fr->step - ener_old->step_prev;
 896         for (i = 0; i < fr->nre; i++)
 897         {
 898             esum_all         = fr->ener[i].esum;
 899             eav_all          = fr->ener[i].eav;
 900             fr->ener[i].esum = esum_all - ener_old->ener_prev[i].esum;
 901             fr->ener[i].eav  = eav_all  - ener_old->ener_prev[i].eav
 902                 - dsqr(ener_old->ener_prev[i].esum/(nstep_all - fr->nsum)
 903                        - esum_all/nstep_all)*
 904                 (nstep_all - fr->nsum)*nstep_all/(double)fr->nsum;
 905             ener_old->ener_prev[i].esum = esum_all;
 906             ener_old->ener_prev[i].eav  = eav_all;
 907         }
 908         ener_old->nsum_prev = nstep_all;
 909     }
 910     else if (fr->nsum > 0)
 911     {
 912         if (fr->nsum != nstep_all)
 913         {
 914             fprintf(stderr, "\nWARNING: something is wrong with the energy sums, will not use exact averages\n");
 915             ener_old->nsum_prev = 0;
 916         }
 917         else
 918         {
 919             ener_old->nsum_prev = nstep_all;
 920         }
 921         /* Copy all sums to ener_prev */
 922         for (i = 0; i < fr->nre; i++)
 923         {
 924             ener_old->ener_prev[i].esum = fr->ener[i].esum;
 925             ener_old->ener_prev[i].eav  = fr->ener[i].eav;
 926         }
 927     }
 928
 929     ener_old->step_prev = fr->step;
 930 }
 931
 932 gmx_bool do_enx(ener_file_t ef, t_enxframe *fr)
 933 {
 934     int           file_version = -1;
 935     int           i, b;
 936     gmx_bool      bRead, bOK, bOK1, bSane;
 937     real          tmp1, tmp2, rdum;
 938     /*int       d_size;*/
 939
 940     bOK   = TRUE;
 941     bRead = gmx_fio_getread(ef->fio);
 942     if (!bRead)
 943     {
 944         fr->e_size = fr->nre*sizeof(fr->ener[0].e)*4;
 945         /*d_size = fr->ndisre*(sizeof(real)*2);*/
 946     }
 947
 948     if (!do_eheader(ef, &file_version, fr, -1, NULL, &bOK))
 949     {
 950         if (bRead)
 951         {
 952             fprintf(stderr, "\rLast energy frame read %d time %8.3f         ",
 953                     ef->framenr-1, ef->frametime);
 954             if (!bOK)
 955             {
 956                 fprintf(stderr,
 957                         "\nWARNING: Incomplete energy frame: nr %d time %8.3f\n",
 958                         ef->framenr, fr->t);
 959             }
 960         }
 961         else
 962         {
 963             gmx_file("Cannot write energy file header; maybe you are out of disk space?");
 964         }
 965         return FALSE;
 966     }
 967     if (bRead)
 968     {
 969         if ((ef->framenr <   20 || ef->framenr %   10 == 0) &&
 970             (ef->framenr <  200 || ef->framenr %  100 == 0) &&
 971             (ef->framenr < 2000 || ef->framenr % 1000 == 0))
 972         {
 973             fprintf(stderr, "\rReading energy frame %6d time %8.3f         ",
 974                     ef->framenr, fr->t);
 975         }
 976         ef->framenr++;
 977         ef->frametime = fr->t;
 978     }
 979     /* Check sanity of this header */
 980     bSane = fr->nre > 0;
 981     for (b = 0; b < fr->nblock; b++)
 982     {
 983         bSane = bSane || (fr->block[b].nsub > 0);
 984     }
 985     if (!((fr->step >= 0) && bSane))
 986     {
 987         fprintf(stderr, "\nWARNING: there may be something wrong with energy file %s\n",
 988                 gmx_fio_getname(ef->fio));
 989         fprintf(stderr, "Found: step=%" GMX_PRId64 ", nre=%d, nblock=%d, time=%g.\n"
 990                 "Trying to skip frame expect a crash though\n",
 991                 fr->step, fr->nre, fr->nblock, fr->t);
 992     }
 993     if (bRead && fr->nre > fr->e_alloc)
 994     {
 995         srenew(fr->ener, fr->nre);
 996         for (i = fr->e_alloc; (i < fr->nre); i++)
 997         {
 998             fr->ener[i].e    = 0;
 999             fr->ener[i].eav  = 0;
1000             fr->ener[i].esum = 0;
1001         }
1002         fr->e_alloc = fr->nre;
1003     }
1004
1005     for (i = 0; i < fr->nre; i++)
1006     {
1007         bOK = bOK && gmx_fio_do_real(ef->fio, fr->ener[i].e);
1008
1009         /* Do not store sums of length 1,
1010          * since this does not add information.
1011          */
1012         if (file_version == 1 ||
1013             (bRead && fr->nsum > 0) || fr->nsum > 1)
1014         {
1015             tmp1 = fr->ener[i].eav;
1016             bOK  = bOK && gmx_fio_do_real(ef->fio, tmp1);
1017             if (bRead)
1018             {
1019                 fr->ener[i].eav = tmp1;
1020             }
1021
1022             /* This is to save only in single precision (unless compiled in DP) */
1023             tmp2 = fr->ener[i].esum;
1024             bOK  = bOK && gmx_fio_do_real(ef->fio, tmp2);
1025             if (bRead)
1026             {
1027                 fr->ener[i].esum = tmp2;
1028             }
1029
1030             if (file_version == 1)
1031             {
1032                 /* Old, unused real */
1033                 rdum = 0;
1034                 bOK  = bOK && gmx_fio_do_real(ef->fio, rdum);
1035             }
1036         }
1037     }
1038
1039     /* Here we can not check for file_version==1, since one could have
1040      * continued an old format simulation with a new one with mdrun -append.
1041      */
1042     if (bRead && ef->eo.bOldFileOpen)
1043     {
1044         /* Convert old full simulation sums to sums between energy frames */
1045         convert_full_sums(&(ef->eo), fr);
1046     }
1047     /* read the blocks */
1048     for (b = 0; b < fr->nblock; b++)
1049     {
1050         /* now read the subblocks. */
1051         int nsub = fr->block[b].nsub; /* shortcut */
1052         int i;
1053
1054         for (i = 0; i < nsub; i++)
1055         {
1056             t_enxsubblock *sub = &(fr->block[b].sub[i]); /* shortcut */
1057
1058             if (bRead)
1059             {
1060                 enxsubblock_alloc(sub);
1061             }
1062
1063             /* read/write data */
1064             switch (sub->type)
1065             {
1066                 case xdr_datatype_float:
1067                     bOK1 = gmx_fio_ndo_float(ef->fio, sub->fval, sub->nr);
1068                     break;
1069                 case xdr_datatype_double:
1070                     bOK1 = gmx_fio_ndo_double(ef->fio, sub->dval, sub->nr);
1071                     break;
1072                 case xdr_datatype_int:
1073                     bOK1 = gmx_fio_ndo_int(ef->fio, sub->ival, sub->nr);
1074                     break;
1075                 case xdr_datatype_int64:
1076                     bOK1 = gmx_fio_ndo_int64(ef->fio, sub->lval, sub->nr);
1077                     break;
1078                 case xdr_datatype_char:
1079                     bOK1 = gmx_fio_ndo_uchar(ef->fio, sub->cval, sub->nr);
1080                     break;
1081                 case xdr_datatype_string:
1082                     bOK1 = gmx_fio_ndo_string(ef->fio, sub->sval, sub->nr);
1083                     break;
1084                 default:
1085                     gmx_incons("Reading unknown block data type: this file is corrupted or from the future");
1086             }
1087             bOK = bOK && bOK1;
1088         }
1089     }
1090
1091     if (!bRead)
1092     {
1093         if (gmx_fio_flush(ef->fio) != 0)
1094         {
1095             gmx_file("Cannot write energy file; maybe you are out of disk space?");
1096         }
1097     }
1098
1099     if (!bOK)
1100     {
1101         if (bRead)
1102         {
1103             fprintf(stderr, "\nLast energy frame read %d",
1104                     ef->framenr-1);
1105             fprintf(stderr, "\nWARNING: Incomplete energy frame: nr %d time %8.3f\n",
1106                     ef->framenr, fr->t);
1107         }
1108         else
1109         {
1110             gmx_fatal(FARGS, "could not write energies");
1111         }
1112         return FALSE;
1113     }
1114
1115     return TRUE;
1116 }
1117
1118 static real find_energy(const char *name, int nre, gmx_enxnm_t *enm,
1119                         t_enxframe *fr)
1120 {
1121     int i;
1122
1123     for (i = 0; i < nre; i++)
1124     {
1125         if (std::strcmp(enm[i].name, name) == 0)
1126         {
1127             return fr->ener[i].e;
1128         }
1129     }
1130
1131     gmx_fatal(FARGS, "Could not find energy term named '%s'", name);
1132
1133     return 0;
1134 }
1135
1136
1137 void get_enx_state(const char *fn, real t, gmx_groups_t *groups, t_inputrec *ir,
1138                    t_state *state)
1139 {
1140     /* Should match the names in mdebin.c */
1141     static const char *boxvel_nm[] = {
1142         "Box-Vel-XX", "Box-Vel-YY", "Box-Vel-ZZ",
1143         "Box-Vel-YX", "Box-Vel-ZX", "Box-Vel-ZY"
1144     };
1145
1146     static const char *baro_nm[] = {
1147         "Barostat"
1148     };
1149
1150
1151     int          ind0[] = { XX, YY, ZZ, YY, ZZ, ZZ };
1152     int          ind1[] = { XX, YY, ZZ, XX, XX, YY };
1153     int          nre, nfr, i, j, ni, npcoupl;
1154     char         buf[STRLEN];
1155     const char  *bufi;
1156     gmx_enxnm_t *enm = NULL;
1157     t_enxframe  *fr;
1158     ener_file_t  in;
1159
1160     in = open_enx(fn, "r");
1161     do_enxnms(in, &nre, &enm);
1162     snew(fr, 1);
1163     nfr = 0;
1164     while ((nfr == 0 || fr->t != t) && do_enx(in, fr))
1165     {
1166         nfr++;
1167     }
1168     close_enx(in);
1169     fprintf(stderr, "\n");
1170
1171     if (nfr == 0 || fr->t != t)
1172     {
1173         gmx_fatal(FARGS, "Could not find frame with time %f in '%s'", t, fn);
1174     }
1175
1176     npcoupl = TRICLINIC(ir->compress) ? 6 : 3;
1177     if (ir->epc == epcPARRINELLORAHMAN)
1178     {
1179         clear_mat(state->boxv);
1180         for (i = 0; i < npcoupl; i++)
1181         {
1182             state->boxv[ind0[i]][ind1[i]] =
1183                 find_energy(boxvel_nm[i], nre, enm, fr);
1184         }
1185         fprintf(stderr, "\nREAD %d BOX VELOCITIES FROM %s\n\n", npcoupl, fn);
1186     }
1187
1188     if (ir->etc == etcNOSEHOOVER)
1189     {
1190         char cns[20];
1191
1192         cns[0] = '\0';
1193
1194         for (i = 0; i < state->ngtc; i++)
1195         {
1196             ni   = groups->grps[egcTC].nm_ind[i];
1197             bufi = *(groups->grpname[ni]);
1198             for (j = 0; (j < state->nhchainlength); j++)
1199             {
1200                 if (IR_NVT_TROTTER(ir))
1201                 {
1202                     sprintf(cns, "-%d", j);
1203                 }
1204                 sprintf(buf, "Xi%s-%s", cns, bufi);
1205                 state->nosehoover_xi[i] = find_energy(buf, nre, enm, fr);
1206                 sprintf(buf, "vXi%s-%s", cns, bufi);
1207                 state->nosehoover_vxi[i] = find_energy(buf, nre, enm, fr);
1208             }
1209
1210         }
1211         fprintf(stderr, "\nREAD %d NOSE-HOOVER Xi chains FROM %s\n\n", state->ngtc, fn);
1212
1213         if (IR_NPT_TROTTER(ir) || IR_NPH_TROTTER(ir))
1214         {
1215             for (i = 0; i < state->nnhpres; i++)
1216             {
1217                 bufi = baro_nm[0]; /* All barostat DOF's together for now */
1218                 for (j = 0; (j < state->nhchainlength); j++)
1219                 {
1220                     sprintf(buf, "Xi-%d-%s", j, bufi);
1221                     state->nhpres_xi[i] = find_energy(buf, nre, enm, fr);
1222                     sprintf(buf, "vXi-%d-%s", j, bufi);
1223                     state->nhpres_vxi[i] = find_energy(buf, nre, enm, fr);
1224                 }
1225             }
1226             fprintf(stderr, "\nREAD %d NOSE-HOOVER BAROSTAT Xi chains FROM %s\n\n", state->nnhpres, fn);
1227         }
1228     }
1229
1230     free_enxnms(nre, enm);
1231     free_enxframe(fr);
1232     sfree(fr);
1233 }