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,2016,2017, 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/math/functions.h"
  50 #include "gromacs/math/vec.h"
  51 #include "gromacs/mdtypes/inputrec.h"
  52 #include "gromacs/mdtypes/md_enums.h"
  53 #include "gromacs/pbcutil/pbc.h"
  54 #include "gromacs/topology/topology.h"
  55 #include "gromacs/utility/compare.h"
  56 #include "gromacs/utility/fatalerror.h"
  57 #include "gromacs/utility/futil.h"
  58 #include "gromacs/utility/gmxassert.h"
  59 #include "gromacs/utility/smalloc.h"
  60
  61 /* The source code in this file should be thread-safe.
  62          Please keep it that way. */
  63
  64 /* This number should be increased whenever the file format changes! */
  65 static const int enx_version = 5;
  66
  67 const char      *enx_block_id_name[] = {
  68     "Averaged orientation restraints",
  69     "Instantaneous orientation restraints",
  70     "Orientation restraint order tensor(s)",
  71     "Distance restraints",
  72     "Free energy data",
  73     "BAR histogram",
  74     "Delta H raw data"
  75 };
  76
  77
  78 /* Stuff for reading pre 4.1 energy files */
  79 typedef struct {
  80     gmx_bool     bOldFileOpen;   /* Is this an open old file? */
  81     gmx_bool     bReadFirstStep; /* Did we read the first step? */
  82     int          first_step;     /* First step in the energy file */
  83     int          step_prev;      /* Previous step */
  84     int          nsum_prev;      /* Previous step sum length */
  85     t_energy    *ener_prev;      /* Previous energy sums */
  86 } ener_old_t;
  87
  88 struct ener_file
  89 {
  90     ener_old_t eo;
  91     t_fileio  *fio;
  92     int        framenr;
  93     real       frametime;
  94 };
  95
  96 static void enxsubblock_init(t_enxsubblock *sb)
  97 {
  98     sb->nr = 0;
  99 #if GMX_DOUBLE
 100     sb->type = xdr_datatype_double;
 101 #else
 102     sb->type = xdr_datatype_float;
 103 #endif
 104     sb->fval       = nullptr;
 105     sb->dval       = nullptr;
 106     sb->ival       = nullptr;
 107     sb->lval       = nullptr;
 108     sb->cval       = nullptr;
 109     sb->sval       = nullptr;
 110     sb->fval_alloc = 0;
 111     sb->dval_alloc = 0;
 112     sb->ival_alloc = 0;
 113     sb->lval_alloc = 0;
 114     sb->cval_alloc = 0;
 115     sb->sval_alloc = 0;
 116 }
 117
 118 static void enxsubblock_free(t_enxsubblock *sb)
 119 {
 120     if (sb->fval_alloc)
 121     {
 122         sfree(sb->fval);
 123         sb->fval_alloc = 0;
 124         sb->fval       = nullptr;
 125     }
 126     if (sb->dval_alloc)
 127     {
 128         sfree(sb->dval);
 129         sb->dval_alloc = 0;
 130         sb->dval       = nullptr;
 131     }
 132     if (sb->ival_alloc)
 133     {
 134         sfree(sb->ival);
 135         sb->ival_alloc = 0;
 136         sb->ival       = nullptr;
 137     }
 138     if (sb->lval_alloc)
 139     {
 140         sfree(sb->lval);
 141         sb->lval_alloc = 0;
 142         sb->lval       = nullptr;
 143     }
 144     if (sb->cval_alloc)
 145     {
 146         sfree(sb->cval);
 147         sb->cval_alloc = 0;
 148         sb->cval       = nullptr;
 149     }
 150     if (sb->sval_alloc)
 151     {
 152         int i;
 153
 154         for (i = 0; i < sb->sval_alloc; i++)
 155         {
 156             if (sb->sval[i])
 157             {
 158                 sfree(sb->sval[i]);
 159             }
 160         }
 161         sfree(sb->sval);
 162         sb->sval_alloc = 0;
 163         sb->sval       = nullptr;
 164     }
 165 }
 166
 167 /* allocate the appropriate amount of memory for the given type and nr */
 168 static void enxsubblock_alloc(t_enxsubblock *sb)
 169 {
 170     /* allocate the appropriate amount of memory */
 171     switch (sb->type)
 172     {
 173         case xdr_datatype_float:
 174             if (sb->nr > sb->fval_alloc)
 175             {
 176                 srenew(sb->fval, sb->nr);
 177                 sb->fval_alloc = sb->nr;
 178             }
 179             break;
 180         case xdr_datatype_double:
 181             if (sb->nr > sb->dval_alloc)
 182             {
 183                 srenew(sb->dval, sb->nr);
 184                 sb->dval_alloc = sb->nr;
 185             }
 186             break;
 187         case xdr_datatype_int:
 188             if (sb->nr > sb->ival_alloc)
 189             {
 190                 srenew(sb->ival, sb->nr);
 191                 sb->ival_alloc = sb->nr;
 192             }
 193             break;
 194         case xdr_datatype_int64:
 195             if (sb->nr > sb->lval_alloc)
 196             {
 197                 srenew(sb->lval, sb->nr);
 198                 sb->lval_alloc = sb->nr;
 199             }
 200             break;
 201         case xdr_datatype_char:
 202             if (sb->nr > sb->cval_alloc)
 203             {
 204                 srenew(sb->cval, sb->nr);
 205                 sb->cval_alloc = sb->nr;
 206             }
 207             break;
 208         case xdr_datatype_string:
 209             if (sb->nr > sb->sval_alloc)
 210             {
 211                 int i;
 212
 213                 srenew(sb->sval, sb->nr);
 214                 for (i = sb->sval_alloc; i < sb->nr; i++)
 215                 {
 216                     sb->sval[i] = nullptr;
 217                 }
 218                 sb->sval_alloc = sb->nr;
 219             }
 220             break;
 221         default:
 222             gmx_incons("Unknown block type: this file is corrupted or from the future");
 223     }
 224 }
 225
 226 static void enxblock_init(t_enxblock *eb)
 227 {
 228     eb->id         = enxOR;
 229     eb->nsub       = 0;
 230     eb->sub        = nullptr;
 231     eb->nsub_alloc = 0;
 232 }
 233
 234 static void enxblock_free(t_enxblock *eb)
 235 {
 236     if (eb->nsub_alloc > 0)
 237     {
 238         int i;
 239         for (i = 0; i < eb->nsub_alloc; i++)
 240         {
 241             enxsubblock_free(&(eb->sub[i]));
 242         }
 243         sfree(eb->sub);
 244         eb->nsub_alloc = 0;
 245         eb->sub        = nullptr;
 246     }
 247 }
 248
 249 void init_enxframe(t_enxframe *fr)
 250 {
 251     fr->e_alloc = 0;
 252     fr->ener    = nullptr;
 253
 254     /*fr->d_alloc=0;*/
 255     /*fr->ndisre=0;*/
 256
 257     fr->nblock       = 0;
 258     fr->nblock_alloc = 0;
 259     fr->block        = nullptr;
 260 }
 261
 262
 263 void free_enxframe(t_enxframe *fr)
 264 {
 265     int b;
 266
 267     if (fr->e_alloc)
 268     {
 269         sfree(fr->ener);
 270     }
 271     for (b = 0; b < fr->nblock_alloc; b++)
 272     {
 273         enxblock_free(&(fr->block[b]));
 274     }
 275     sfree(fr->block);
 276 }
 277
 278 void add_blocks_enxframe(t_enxframe *fr, int n)
 279 {
 280     fr->nblock = n;
 281     if (n > fr->nblock_alloc)
 282     {
 283         int b;
 284
 285         srenew(fr->block, n);
 286         for (b = fr->nblock_alloc; b < fr->nblock; b++)
 287         {
 288             enxblock_init(&(fr->block[b]));
 289         }
 290         fr->nblock_alloc = n;
 291     }
 292 }
 293
 294 t_enxblock *find_block_id_enxframe(t_enxframe *ef, int id, t_enxblock *prev)
 295 {
 296     gmx_off_t starti = 0;
 297     gmx_off_t i;
 298
 299     if (prev)
 300     {
 301         starti = (prev - ef->block) + 1;
 302     }
 303     for (i = starti; i < ef->nblock; i++)
 304     {
 305         if (ef->block[i].id == id)
 306         {
 307             return &(ef->block[i]);
 308         }
 309     }
 310     return nullptr;
 311 }
 312
 313 void add_subblocks_enxblock(t_enxblock *eb, int n)
 314 {
 315     eb->nsub = n;
 316     if (eb->nsub > eb->nsub_alloc)
 317     {
 318         int b;
 319
 320         srenew(eb->sub, n);
 321         for (b = eb->nsub_alloc; b < n; b++)
 322         {
 323             enxsubblock_init(&(eb->sub[b]));
 324         }
 325         eb->nsub_alloc = n;
 326     }
 327 }
 328
 329 static void enx_warning(const char *msg)
 330 {
 331     if (getenv("GMX_ENX_NO_FATAL") != nullptr)
 332     {
 333         gmx_warning(msg);
 334     }
 335     else
 336     {
 337         gmx_fatal(FARGS, "%s\n%s",
 338                   msg,
 339                   "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");
 340     }
 341 }
 342
 343 static void edr_strings(XDR *xdr, gmx_bool bRead, int file_version,
 344                         int n, gmx_enxnm_t **nms)
 345 {
 346     int          i;
 347     gmx_enxnm_t *nm;
 348
 349     if (*nms == nullptr)
 350     {
 351         snew(*nms, n);
 352     }
 353     for (i = 0; i < n; i++)
 354     {
 355         nm = &(*nms)[i];
 356         if (bRead)
 357         {
 358             if (nm->name)
 359             {
 360                 sfree(nm->name);
 361                 nm->name = nullptr;
 362             }
 363             if (nm->unit)
 364             {
 365                 sfree(nm->unit);
 366                 nm->unit = nullptr;
 367             }
 368         }
 369         if (!xdr_string(xdr, &(nm->name), STRLEN))
 370         {
 371             gmx_file("Cannot write energy names to file; maybe you are out of disk space?");
 372         }
 373         if (file_version >= 2)
 374         {
 375             if (!xdr_string(xdr, &(nm->unit), STRLEN))
 376             {
 377                 gmx_file("Cannot write energy names to file; maybe you are out of disk space?");
 378             }
 379         }
 380         else
 381         {
 382             nm->unit = gmx_strdup("kJ/mol");
 383         }
 384     }
 385 }
 386
 387 void do_enxnms(ener_file_t ef, int *nre, gmx_enxnm_t **nms)
 388 {
 389     int      magic = -55555;
 390     XDR     *xdr;
 391     gmx_bool bRead = gmx_fio_getread(ef->fio);
 392     int      file_version;
 393
 394     xdr = gmx_fio_getxdr(ef->fio);
 395
 396     if (!xdr_int(xdr, &magic))
 397     {
 398         if (!bRead)
 399         {
 400             gmx_file("Cannot write energy names to file; maybe you are out of disk space?");
 401         }
 402         *nre = 0;
 403         return;
 404     }
 405     if (magic > 0)
 406     {
 407         /* Assume this is an old edr format */
 408         file_version          = 1;
 409         *nre                  = magic;
 410         ef->eo.bOldFileOpen   = TRUE;
 411         ef->eo.bReadFirstStep = FALSE;
 412         srenew(ef->eo.ener_prev, *nre);
 413     }
 414     else
 415     {
 416         ef->eo.bOldFileOpen = FALSE;
 417
 418         if (magic != -55555)
 419         {
 420             gmx_fatal(FARGS, "Energy names magic number mismatch, this is not a GROMACS edr file");
 421         }
 422         file_version = enx_version;
 423         xdr_int(xdr, &file_version);
 424         if (file_version > enx_version)
 425         {
 426             gmx_fatal(FARGS, "reading tpx file (%s) version %d with version %d program", gmx_fio_getname(ef->fio), file_version, enx_version);
 427         }
 428         xdr_int(xdr, nre);
 429     }
 430     if (file_version != enx_version)
 431     {
 432         fprintf(stderr, "Note: enx file_version %d, software version %d\n",
 433                 file_version, enx_version);
 434     }
 435
 436     edr_strings(xdr, bRead, file_version, *nre, nms);
 437 }
 438
 439 static gmx_bool do_eheader(ener_file_t ef, int *file_version, t_enxframe *fr,
 440                            int nre_test, gmx_bool *bWrongPrecision, gmx_bool *bOK)
 441 {
 442     int          magic = -7777777;
 443     real         first_real_to_check;
 444     int          b, zero = 0, dum = 0;
 445     gmx_bool     bRead      = gmx_fio_getread(ef->fio);
 446     int          ndisre     = 0;
 447     int          startb     = 0;
 448 #if !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         // cppcheck-suppress redundantPointerOp
 500         if (!gmx_fio_do_int(ef->fio, *file_version))
 501         {
 502             *bOK = FALSE;
 503         }
 504         if (*bOK && *file_version > enx_version)
 505         {
 506             gmx_fatal(FARGS, "reading tpx file (%s) version %d with version %d program", gmx_fio_getname(ef->fio), file_version, enx_version);
 507         }
 508         if (!gmx_fio_do_double(ef->fio, fr->t))
 509         {
 510             *bOK = FALSE;
 511         }
 512         if (!gmx_fio_do_int64(ef->fio, fr->step))
 513         {
 514             *bOK = FALSE;
 515         }
 516         if (!bRead && fr->nsum == 1)
 517         {
 518             /* Do not store sums of length 1,
 519              * since this does not add information.
 520              */
 521             if (!gmx_fio_do_int(ef->fio, zero))
 522             {
 523                 *bOK = FALSE;
 524             }
 525         }
 526         else
 527         {
 528             if (!gmx_fio_do_int(ef->fio, fr->nsum))
 529             {
 530                 *bOK = FALSE;
 531             }
 532         }
 533         if (*file_version >= 3)
 534         {
 535             if (!gmx_fio_do_int64(ef->fio, fr->nsteps))
 536             {
 537                 *bOK = FALSE;
 538             }
 539         }
 540         else
 541         {
 542             fr->nsteps = std::max(1, fr->nsum);
 543         }
 544         if (*file_version >= 5)
 545         {
 546             if (!gmx_fio_do_double(ef->fio, fr->dt))
 547             {
 548                 *bOK = FALSE;
 549             }
 550         }
 551         else
 552         {
 553             fr->dt = 0;
 554         }
 555     }
 556     if (!gmx_fio_do_int(ef->fio, fr->nre))
 557     {
 558         *bOK = FALSE;
 559     }
 560     if (*file_version < 4)
 561     {
 562         if (!gmx_fio_do_int(ef->fio, ndisre))
 563         {
 564             *bOK = FALSE;
 565         }
 566     }
 567     else
 568     {
 569         /* now reserved for possible future use */
 570         if (!gmx_fio_do_int(ef->fio, dum))
 571         {
 572             *bOK = FALSE;
 573         }
 574     }
 575
 576     if (!gmx_fio_do_int(ef->fio, fr->nblock))
 577     {
 578         *bOK = FALSE;
 579     }
 580     if (fr->nblock < 0)
 581     {
 582         *bOK = FALSE;
 583     }
 584
 585     if (ndisre != 0)
 586     {
 587         if (*file_version >= 4)
 588         {
 589             enx_warning("Distance restraint blocks in old style in new style file");
 590             *bOK = FALSE;
 591             return FALSE;
 592         }
 593         fr->nblock += 1;
 594     }
 595
 596
 597     /* Frames could have nre=0, so we can not rely only on the fr->nre check */
 598     if (bRead && nre_test >= 0 &&
 599         ((fr->nre > 0 && fr->nre != nre_test) ||
 600          fr->nre < 0 || ndisre < 0 || fr->nblock < 0))
 601     {
 602         if (bWrongPrecision)
 603         {
 604             *bWrongPrecision = TRUE;
 605         }
 606         return *bOK;
 607     }
 608
 609     /* we now know what these should be, or we've already bailed out because
 610        of wrong precision */
 611     if (*file_version == 1 && (fr->t < 0 || fr->t > 1e20 || fr->step < 0 ) )
 612     {
 613         enx_warning("edr file with negative step number or unreasonable time (and without version number).");
 614         *bOK = FALSE;
 615         return FALSE;
 616     }
 617
 618
 619     if (*bOK && bRead)
 620     {
 621         add_blocks_enxframe(fr, fr->nblock);
 622     }
 623
 624     startb = 0;
 625     if (ndisre > 0)
 626     {
 627         /* sub[0] is the instantaneous data, sub[1] is time averaged */
 628         add_subblocks_enxblock(&(fr->block[0]), 2);
 629         fr->block[0].id          = enxDISRE;
 630         fr->block[0].sub[0].nr   = ndisre;
 631         fr->block[0].sub[1].nr   = ndisre;
 632         fr->block[0].sub[0].type = dtreal;
 633         fr->block[0].sub[1].type = dtreal;
 634         startb++;
 635     }
 636
 637     /* read block header info */
 638     for (b = startb; b < fr->nblock; b++)
 639     {
 640         if (*file_version < 4)
 641         {
 642             /* blocks in old version files always have 1 subblock that
 643                consists of reals. */
 644             int nrint;
 645
 646             if (bRead)
 647             {
 648                 add_subblocks_enxblock(&(fr->block[b]), 1);
 649             }
 650             else
 651             {
 652                 if (fr->block[b].nsub != 1)
 653                 {
 654                     gmx_incons("Writing an old version .edr file with too many subblocks");
 655                 }
 656                 if (fr->block[b].sub[0].type != dtreal)
 657                 {
 658                     gmx_incons("Writing an old version .edr file the wrong subblock type");
 659                 }
 660             }
 661             nrint = fr->block[b].sub[0].nr;
 662
 663             if (!gmx_fio_do_int(ef->fio, nrint))
 664             {
 665                 *bOK = FALSE;
 666             }
 667             fr->block[b].id          = b - startb;
 668             fr->block[b].sub[0].nr   = nrint;
 669             fr->block[b].sub[0].type = dtreal;
 670         }
 671         else
 672         {
 673             int i;
 674             /* in the new version files, the block header only contains
 675                the ID and the number of subblocks */
 676             int nsub = fr->block[b].nsub;
 677             *bOK = *bOK && gmx_fio_do_int(ef->fio, fr->block[b].id);
 678             *bOK = *bOK && gmx_fio_do_int(ef->fio, nsub);
 679
 680             fr->block[b].nsub = nsub;
 681             if (bRead)
 682             {
 683                 add_subblocks_enxblock(&(fr->block[b]), nsub);
 684             }
 685
 686             /* read/write type & size for each subblock */
 687             for (i = 0; i < nsub; i++)
 688             {
 689                 t_enxsubblock *sub    = &(fr->block[b].sub[i]); /* shortcut */
 690                 int            typenr = sub->type;
 691
 692                 *bOK = *bOK && gmx_fio_do_int(ef->fio, typenr);
 693                 *bOK = *bOK && gmx_fio_do_int(ef->fio, sub->nr);
 694
 695                 sub->type = (xdr_datatype)typenr;
 696             }
 697         }
 698     }
 699     if (!gmx_fio_do_int(ef->fio, fr->e_size))
 700     {
 701         *bOK = FALSE;
 702     }
 703
 704     /* now reserved for possible future use */
 705     if (!gmx_fio_do_int(ef->fio, dum))
 706     {
 707         *bOK = FALSE;
 708     }
 709
 710     /* Do a dummy int to keep the format compatible with the old code */
 711     if (!gmx_fio_do_int(ef->fio, dum))
 712     {
 713         *bOK = FALSE;
 714     }
 715
 716     if (*bOK && *file_version == 1 && nre_test < 0)
 717     {
 718         if (!ef->eo.bReadFirstStep)
 719         {
 720             ef->eo.bReadFirstStep = TRUE;
 721             ef->eo.first_step     = fr->step;
 722             ef->eo.step_prev      = fr->step;
 723             ef->eo.nsum_prev      = 0;
 724         }
 725
 726         fr->nsum   = fr->step - ef->eo.first_step + 1;
 727         fr->nsteps = fr->step - ef->eo.step_prev;
 728         fr->dt     = 0;
 729     }
 730
 731     return *bOK;
 732 }
 733
 734 void free_enxnms(int n, gmx_enxnm_t *nms)
 735 {
 736     int i;
 737
 738     for (i = 0; i < n; i++)
 739     {
 740         sfree(nms[i].name);
 741         sfree(nms[i].unit);
 742     }
 743
 744     sfree(nms);
 745 }
 746
 747 void close_enx(ener_file_t ef)
 748 {
 749     if (ef == nullptr)
 750     {
 751         // Nothing to do
 752         return;
 753     }
 754     if (gmx_fio_close(ef->fio) != 0)
 755     {
 756         gmx_file("Cannot close energy file; it might be corrupt, or maybe you are out of disk space?");
 757     }
 758 }
 759
 760 void done_ener_file(ener_file_t ef)
 761 {
 762     // Free the contents, then the pointer itself
 763     close_enx(ef);
 764     sfree(ef);
 765 }
 766
 767 /*!\brief Return TRUE if a file exists but is empty, otherwise FALSE.
 768  *
 769  * If the file exists but has length larger than zero, if it does not exist, or
 770  * if there is a file system error, FALSE will be returned instead.
 771  *
 772  * \param fn  File name to test
 773  *
 774  * \return TRUE if file could be open but is empty, otherwise FALSE.
 775  */
 776 static gmx_bool
 777 empty_file(const char *fn)
 778 {
 779     FILE    *fp;
 780     char     dum;
 781     int      ret;
 782     gmx_bool bEmpty;
 783
 784     fp     = gmx_fio_fopen(fn, "r");
 785     ret    = fread(&dum, sizeof(dum), 1, fp);
 786     bEmpty = feof(fp);
 787     gmx_fio_fclose(fp);
 788
 789     // bEmpty==TRUE but ret!=0 would likely be some strange I/O error, but at
 790     // least it's not a normal empty file, so we return FALSE in that case.
 791     return (bEmpty && ret == 0);
 792 }
 793
 794
 795 ener_file_t open_enx(const char *fn, const char *mode)
 796 {
 797     int               nre;
 798     gmx_enxnm_t      *nms          = nullptr;
 799     int               file_version = -1;
 800     t_enxframe       *fr;
 801     gmx_bool          bWrongPrecision, bOK = TRUE;
 802     struct ener_file *ef;
 803
 804     snew(ef, 1);
 805
 806     if (mode[0] == 'r')
 807     {
 808         ef->fio = gmx_fio_open(fn, mode);
 809         gmx_fio_setprecision(ef->fio, FALSE);
 810         do_enxnms(ef, &nre, &nms);
 811         snew(fr, 1);
 812         do_eheader(ef, &file_version, fr, nre, &bWrongPrecision, &bOK);
 813         if (!bOK)
 814         {
 815             gmx_file("Cannot read energy file header. Corrupt file?");
 816         }
 817
 818         /* Now check whether this file is in single precision */
 819         if (!bWrongPrecision &&
 820             ((fr->e_size && (fr->nre == nre) &&
 821               (nre*4*(long int)sizeof(float) == fr->e_size)) ) )
 822         {
 823             fprintf(stderr, "Opened %s as single precision energy file\n", fn);
 824             free_enxnms(nre, nms);
 825         }
 826         else
 827         {
 828             gmx_fio_rewind(ef->fio);
 829             gmx_fio_setprecision(ef->fio, TRUE);
 830             do_enxnms(ef, &nre, &nms);
 831             do_eheader(ef, &file_version, fr, nre, &bWrongPrecision, &bOK);
 832             if (!bOK)
 833             {
 834                 gmx_file("Cannot write energy file header; maybe you are out of disk space?");
 835             }
 836
 837             if (((fr->e_size && (fr->nre == nre) &&
 838                   (nre*4*(long int)sizeof(double) == fr->e_size)) ))
 839             {
 840                 fprintf(stderr, "Opened %s as double precision energy file\n",
 841                         fn);
 842             }
 843             else
 844             {
 845                 if (empty_file(fn))
 846                 {
 847                     gmx_fatal(FARGS, "File %s is empty", fn);
 848                 }
 849                 else
 850                 {
 851                     gmx_fatal(FARGS, "Energy file %s not recognized, maybe different CPU?",
 852                               fn);
 853                 }
 854             }
 855             free_enxnms(nre, nms);
 856         }
 857         free_enxframe(fr);
 858         sfree(fr);
 859         gmx_fio_rewind(ef->fio);
 860     }
 861     else
 862     {
 863         ef->fio = gmx_fio_open(fn, mode);
 864     }
 865
 866     ef->framenr   = 0;
 867     ef->frametime = 0;
 868     return ef;
 869 }
 870
 871 t_fileio *enx_file_pointer(const ener_file_t ef)
 872 {
 873     return ef->fio;
 874 }
 875
 876 static void convert_full_sums(ener_old_t *ener_old, t_enxframe *fr)
 877 {
 878     int    nstep_all;
 879     int    ne, ns, i;
 880     double esum_all, eav_all;
 881
 882     if (fr->nsum > 0)
 883     {
 884         ne = 0;
 885         ns = 0;
 886         for (i = 0; i < fr->nre; i++)
 887         {
 888             if (fr->ener[i].e    != 0)
 889             {
 890                 ne++;
 891             }
 892             if (fr->ener[i].esum != 0)
 893             {
 894                 ns++;
 895             }
 896         }
 897         if (ne > 0 && ns == 0)
 898         {
 899             /* We do not have all energy sums */
 900             fr->nsum = 0;
 901         }
 902     }
 903
 904     /* Convert old full simulation sums to sums between energy frames */
 905     nstep_all = fr->step - ener_old->first_step + 1;
 906     if (fr->nsum > 1 && fr->nsum == nstep_all && ener_old->nsum_prev > 0)
 907     {
 908         /* Set the new sum length: the frame step difference */
 909         fr->nsum = fr->step - ener_old->step_prev;
 910         for (i = 0; i < fr->nre; i++)
 911         {
 912             esum_all         = fr->ener[i].esum;
 913             eav_all          = fr->ener[i].eav;
 914             fr->ener[i].esum = esum_all - ener_old->ener_prev[i].esum;
 915             fr->ener[i].eav  = eav_all  - ener_old->ener_prev[i].eav
 916                 - gmx::square(ener_old->ener_prev[i].esum/(nstep_all - fr->nsum)
 917                               - esum_all/nstep_all)*
 918                 (nstep_all - fr->nsum)*nstep_all/(double)fr->nsum;
 919             ener_old->ener_prev[i].esum = esum_all;
 920             ener_old->ener_prev[i].eav  = eav_all;
 921         }
 922         ener_old->nsum_prev = nstep_all;
 923     }
 924     else if (fr->nsum > 0)
 925     {
 926         if (fr->nsum != nstep_all)
 927         {
 928             fprintf(stderr, "\nWARNING: something is wrong with the energy sums, will not use exact averages\n");
 929             ener_old->nsum_prev = 0;
 930         }
 931         else
 932         {
 933             ener_old->nsum_prev = nstep_all;
 934         }
 935         /* Copy all sums to ener_prev */
 936         for (i = 0; i < fr->nre; i++)
 937         {
 938             ener_old->ener_prev[i].esum = fr->ener[i].esum;
 939             ener_old->ener_prev[i].eav  = fr->ener[i].eav;
 940         }
 941     }
 942
 943     ener_old->step_prev = fr->step;
 944 }
 945
 946 gmx_bool do_enx(ener_file_t ef, t_enxframe *fr)
 947 {
 948     int           file_version = -1;
 949     int           i, b;
 950     gmx_bool      bRead, bOK, bOK1, bSane;
 951     real          tmp1, tmp2, rdum;
 952     /*int       d_size;*/
 953
 954     bOK   = TRUE;
 955     bRead = gmx_fio_getread(ef->fio);
 956     if (!bRead)
 957     {
 958         fr->e_size = fr->nre*sizeof(fr->ener[0].e)*4;
 959         /*d_size = fr->ndisre*(sizeof(real)*2);*/
 960     }
 961
 962     if (!do_eheader(ef, &file_version, fr, -1, nullptr, &bOK))
 963     {
 964         if (bRead)
 965         {
 966             fprintf(stderr, "\rLast energy frame read %d time %8.3f         ",
 967                     ef->framenr-1, ef->frametime);
 968             fflush(stderr);
 969
 970             if (!bOK)
 971             {
 972                 fprintf(stderr,
 973                         "\nWARNING: Incomplete energy frame: nr %d time %8.3f\n",
 974                         ef->framenr, fr->t);
 975             }
 976         }
 977         else
 978         {
 979             gmx_file("Cannot write energy file header; maybe you are out of disk space?");
 980         }
 981         return FALSE;
 982     }
 983     if (bRead)
 984     {
 985         if ((ef->framenr <   20 || ef->framenr %   10 == 0) &&
 986             (ef->framenr <  200 || ef->framenr %  100 == 0) &&
 987             (ef->framenr < 2000 || ef->framenr % 1000 == 0))
 988         {
 989             fprintf(stderr, "\rReading energy frame %6d time %8.3f         ",
 990                     ef->framenr, fr->t);
 991         }
 992         ef->framenr++;
 993         ef->frametime = fr->t;
 994     }
 995     /* Check sanity of this header */
 996     bSane = fr->nre > 0;
 997     for (b = 0; b < fr->nblock; b++)
 998     {
 999         bSane = bSane || (fr->block[b].nsub > 0);
1000     }
1001     if (!((fr->step >= 0) && bSane) && bRead)
1002     {
1003         fprintf(stderr, "\nWARNING: there may be something wrong with energy file %s\n",
1004                 gmx_fio_getname(ef->fio));
1005         fprintf(stderr, "Found: step=%" GMX_PRId64 ", nre=%d, nblock=%d, time=%g.\n",
1006                 fr->step, fr->nre, fr->nblock, fr->t);
1007     }
1008     if (bRead && fr->nre > fr->e_alloc)
1009     {
1010         srenew(fr->ener, fr->nre);
1011         for (i = fr->e_alloc; (i < fr->nre); i++)
1012         {
1013             fr->ener[i].e    = 0;
1014             fr->ener[i].eav  = 0;
1015             fr->ener[i].esum = 0;
1016         }
1017         fr->e_alloc = fr->nre;
1018     }
1019
1020     for (i = 0; i < fr->nre; i++)
1021     {
1022         bOK = bOK && gmx_fio_do_real(ef->fio, fr->ener[i].e);
1023
1024         /* Do not store sums of length 1,
1025          * since this does not add information.
1026          */
1027         if (file_version == 1 ||
1028             (bRead && fr->nsum > 0) || fr->nsum > 1)
1029         {
1030             tmp1 = fr->ener[i].eav;
1031             bOK  = bOK && gmx_fio_do_real(ef->fio, tmp1);
1032             if (bRead)
1033             {
1034                 fr->ener[i].eav = tmp1;
1035             }
1036
1037             /* This is to save only in single precision (unless compiled in DP) */
1038             tmp2 = fr->ener[i].esum;
1039             bOK  = bOK && gmx_fio_do_real(ef->fio, tmp2);
1040             if (bRead)
1041             {
1042                 fr->ener[i].esum = tmp2;
1043             }
1044
1045             if (file_version == 1)
1046             {
1047                 /* Old, unused real */
1048                 rdum = 0;
1049                 bOK  = bOK && gmx_fio_do_real(ef->fio, rdum);
1050             }
1051         }
1052     }
1053
1054     /* Here we can not check for file_version==1, since one could have
1055      * continued an old format simulation with a new one with mdrun -append.
1056      */
1057     if (bRead && ef->eo.bOldFileOpen)
1058     {
1059         /* Convert old full simulation sums to sums between energy frames */
1060         convert_full_sums(&(ef->eo), fr);
1061     }
1062     /* read the blocks */
1063     for (b = 0; b < fr->nblock; b++)
1064     {
1065         /* now read the subblocks. */
1066         int nsub = fr->block[b].nsub; /* shortcut */
1067         int i;
1068
1069         for (i = 0; i < nsub; i++)
1070         {
1071             t_enxsubblock *sub = &(fr->block[b].sub[i]); /* shortcut */
1072
1073             if (bRead)
1074             {
1075                 enxsubblock_alloc(sub);
1076             }
1077
1078             /* read/write data */
1079             switch (sub->type)
1080             {
1081                 case xdr_datatype_float:
1082                     bOK1 = gmx_fio_ndo_float(ef->fio, sub->fval, sub->nr);
1083                     break;
1084                 case xdr_datatype_double:
1085                     bOK1 = gmx_fio_ndo_double(ef->fio, sub->dval, sub->nr);
1086                     break;
1087                 case xdr_datatype_int:
1088                     bOK1 = gmx_fio_ndo_int(ef->fio, sub->ival, sub->nr);
1089                     break;
1090                 case xdr_datatype_int64:
1091                     bOK1 = gmx_fio_ndo_int64(ef->fio, sub->lval, sub->nr);
1092                     break;
1093                 case xdr_datatype_char:
1094                     bOK1 = gmx_fio_ndo_uchar(ef->fio, sub->cval, sub->nr);
1095                     break;
1096                 case xdr_datatype_string:
1097                     bOK1 = gmx_fio_ndo_string(ef->fio, sub->sval, sub->nr);
1098                     break;
1099                 default:
1100                     gmx_incons("Reading unknown block data type: this file is corrupted or from the future");
1101             }
1102             bOK = bOK && bOK1;
1103         }
1104     }
1105
1106     if (!bRead)
1107     {
1108         if (gmx_fio_flush(ef->fio) != 0)
1109         {
1110             gmx_file("Cannot write energy file; maybe you are out of disk space?");
1111         }
1112     }
1113
1114     if (!bOK)
1115     {
1116         if (bRead)
1117         {
1118             fprintf(stderr, "\nLast energy frame read %d",
1119                     ef->framenr-1);
1120             fprintf(stderr, "\nWARNING: Incomplete energy frame: nr %d time %8.3f\n",
1121                     ef->framenr, fr->t);
1122         }
1123         else
1124         {
1125             gmx_fatal(FARGS, "could not write energies");
1126         }
1127         return FALSE;
1128     }
1129
1130     return TRUE;
1131 }
1132
1133 static real find_energy(const char *name, int nre, gmx_enxnm_t *enm,
1134                         t_enxframe *fr)
1135 {
1136     int i;
1137
1138     for (i = 0; i < nre; i++)
1139     {
1140         if (std::strcmp(enm[i].name, name) == 0)
1141         {
1142             return fr->ener[i].e;
1143         }
1144     }
1145
1146     gmx_fatal(FARGS, "Could not find energy term named '%s'", name);
1147
1148     return 0;
1149 }
1150
1151
1152 void get_enx_state(const char *fn, real t, const gmx_groups_t *groups, t_inputrec *ir,
1153                    t_state *state)
1154 {
1155     /* Should match the names in mdebin.c */
1156     static const char *boxvel_nm[] = {
1157         "Box-Vel-XX", "Box-Vel-YY", "Box-Vel-ZZ",
1158         "Box-Vel-YX", "Box-Vel-ZX", "Box-Vel-ZY"
1159     };
1160
1161     static const char *baro_nm[] = {
1162         "Barostat"
1163     };
1164
1165
1166     int          ind0[] = { XX, YY, ZZ, YY, ZZ, ZZ };
1167     int          ind1[] = { XX, YY, ZZ, XX, XX, YY };
1168     int          nre, nfr, i, j, ni, npcoupl;
1169     char         buf[STRLEN];
1170     const char  *bufi;
1171     gmx_enxnm_t *enm = nullptr;
1172     t_enxframe  *fr;
1173     ener_file_t  in;
1174
1175     in = open_enx(fn, "r");
1176     do_enxnms(in, &nre, &enm);
1177     snew(fr, 1);
1178     nfr = 0;
1179     while ((nfr == 0 || fr->t != t) && do_enx(in, fr))
1180     {
1181         nfr++;
1182     }
1183     close_enx(in);
1184     fprintf(stderr, "\n");
1185
1186     if (nfr == 0 || fr->t != t)
1187     {
1188         gmx_fatal(FARGS, "Could not find frame with time %f in '%s'", t, fn);
1189     }
1190
1191     npcoupl = TRICLINIC(ir->compress) ? 6 : 3;
1192     if (ir->epc == epcPARRINELLORAHMAN)
1193     {
1194         clear_mat(state->boxv);
1195         for (i = 0; i < npcoupl; i++)
1196         {
1197             state->boxv[ind0[i]][ind1[i]] =
1198                 find_energy(boxvel_nm[i], nre, enm, fr);
1199         }
1200         fprintf(stderr, "\nREAD %d BOX VELOCITIES FROM %s\n\n", npcoupl, fn);
1201     }
1202
1203     if (ir->etc == etcNOSEHOOVER)
1204     {
1205         char cns[20];
1206
1207         cns[0] = '\0';
1208
1209         for (i = 0; i < state->ngtc; i++)
1210         {
1211             ni   = groups->grps[egcTC].nm_ind[i];
1212             bufi = *(groups->grpname[ni]);
1213             for (j = 0; (j < state->nhchainlength); j++)
1214             {
1215                 if (inputrecNvtTrotter(ir))
1216                 {
1217                     sprintf(cns, "-%d", j);
1218                 }
1219                 sprintf(buf, "Xi%s-%s", cns, bufi);
1220                 state->nosehoover_xi[i] = find_energy(buf, nre, enm, fr);
1221                 sprintf(buf, "vXi%s-%s", cns, bufi);
1222                 state->nosehoover_vxi[i] = find_energy(buf, nre, enm, fr);
1223             }
1224
1225         }
1226         fprintf(stderr, "\nREAD %d NOSE-HOOVER Xi chains FROM %s\n\n", state->ngtc, fn);
1227
1228         if (inputrecNptTrotter(ir) || inputrecNphTrotter(ir))
1229         {
1230             for (i = 0; i < state->nnhpres; i++)
1231             {
1232                 bufi = baro_nm[0]; /* All barostat DOF's together for now */
1233                 for (j = 0; (j < state->nhchainlength); j++)
1234                 {
1235                     sprintf(buf, "Xi-%d-%s", j, bufi);
1236                     state->nhpres_xi[i] = find_energy(buf, nre, enm, fr);
1237                     sprintf(buf, "vXi-%d-%s", j, bufi);
1238                     state->nhpres_vxi[i] = find_energy(buf, nre, enm, fr);
1239                 }
1240             }
1241             fprintf(stderr, "\nREAD %d NOSE-HOOVER BAROSTAT Xi chains FROM %s\n\n", state->nnhpres, fn);
1242         }
1243     }
1244
1245     free_enxnms(nre, enm);
1246     free_enxframe(fr);
1247     sfree(fr);
1248 }
1249
1250 static real ener_tensor_diag(int n, int *ind1, int *ind2,
1251                              gmx_enxnm_t *enm1,
1252                              int *tensi, int i,
1253                              t_energy e1[], t_energy e2[])
1254 {
1255     int    d1, d2;
1256     int    j;
1257     real   prod1, prod2;
1258     int    nfound;
1259     size_t len;
1260
1261     d1 = tensi[i]/DIM;
1262     d2 = tensi[i] - d1*DIM;
1263
1264     /* Find the diagonal elements d1 and d2 */
1265     len    = std::strlen(enm1[ind1[i]].name);
1266     prod1  = 1;
1267     prod2  = 1;
1268     nfound = 0;
1269     for (j = 0; j < n; j++)
1270     {
1271         if (tensi[j] >= 0 &&
1272             std::strlen(enm1[ind1[j]].name) == len &&
1273             std::strncmp(enm1[ind1[i]].name, enm1[ind1[j]].name, len-2) == 0 &&
1274             (tensi[j] == d1*DIM+d1 || tensi[j] == d2*DIM+d2))
1275         {
1276             prod1 *= fabs(e1[ind1[j]].e);
1277             prod2 *= fabs(e2[ind2[j]].e);
1278             nfound++;
1279         }
1280     }
1281
1282     if (nfound == 2)
1283     {
1284         return 0.5*(std::sqrt(prod1) + std::sqrt(prod2));
1285     }
1286     else
1287     {
1288         return 0;
1289     }
1290 }
1291
1292 static gmx_bool enernm_equal(const char *nm1, const char *nm2)
1293 {
1294     int len1, len2;
1295
1296     len1 = std::strlen(nm1);
1297     len2 = std::strlen(nm2);
1298
1299     /* Remove " (bar)" at the end of a name */
1300     if (len1 > 6 && std::strcmp(nm1+len1-6, " (bar)") == 0)
1301     {
1302         len1 -= 6;
1303     }
1304     if (len2 > 6 && std::strcmp(nm2+len2-6, " (bar)") == 0)
1305     {
1306         len2 -= 6;
1307     }
1308
1309     return (len1 == len2 && gmx_strncasecmp(nm1, nm2, len1) == 0);
1310 }
1311
1312 static void cmp_energies(FILE *fp, int step1, int step2,
1313                          t_energy e1[], t_energy e2[],
1314                          gmx_enxnm_t *enm1,
1315                          real ftol, real abstol,
1316                          int nre, int *ind1, int *ind2, int maxener)
1317 {
1318     int   i, ii;
1319     int  *tensi, len, d1, d2;
1320     real  ftol_i, abstol_i;
1321
1322     snew(tensi, maxener);
1323     /* Check for tensor elements ending on "-XX", "-XY", ... , "-ZZ" */
1324     for (i = 0; (i < maxener); i++)
1325     {
1326         ii       = ind1[i];
1327         tensi[i] = -1;
1328         len      = std::strlen(enm1[ii].name);
1329         if (len > 3 && enm1[ii].name[len-3] == '-')
1330         {
1331             d1 = enm1[ii].name[len-2] - 'X';
1332             d2 = enm1[ii].name[len-1] - 'X';
1333             if (d1 >= 0 && d1 < DIM &&
1334                 d2 >= 0 && d2 < DIM)
1335             {
1336                 tensi[i] = d1*DIM + d2;
1337             }
1338         }
1339     }
1340
1341     for (i = 0; (i < maxener); i++)
1342     {
1343         /* Check if this is an off-diagonal tensor element */
1344         if (tensi[i] >= 0 && tensi[i] != 0 && tensi[i] != 4 && tensi[i] != 8)
1345         {
1346             /* Turn on the relative tolerance check (4 is maximum relative diff.) */
1347             ftol_i = 5;
1348             /* Do the relative tolerance through an absolute tolerance times
1349              * the size of diagonal components of the tensor.
1350              */
1351             abstol_i = ftol*ener_tensor_diag(nre, ind1, ind2, enm1, tensi, i, e1, e2);
1352             if (debug)
1353             {
1354                 fprintf(debug, "tensor '%s' val %f diag %f\n",
1355                         enm1[i].name, e1[i].e, abstol_i/ftol);
1356             }
1357             if (abstol_i > 0)
1358             {
1359                 /* We found a diagonal, we need to check with the minimum tolerance */
1360                 abstol_i = std::min(abstol_i, abstol);
1361             }
1362             else
1363             {
1364                 /* We did not find a diagonal, ignore the relative tolerance check */
1365                 abstol_i = abstol;
1366             }
1367         }
1368         else
1369         {
1370             ftol_i   = ftol;
1371             abstol_i = abstol;
1372         }
1373         if (!equal_real(e1[ind1[i]].e, e2[ind2[i]].e, ftol_i, abstol_i))
1374         {
1375             fprintf(fp, "%-15s  step %3d:  %12g,  step %3d: %12g\n",
1376                     enm1[ind1[i]].name,
1377                     step1, e1[ind1[i]].e,
1378                     step2, e2[ind2[i]].e);
1379         }
1380     }
1381
1382     sfree(tensi);
1383 }
1384
1385 #if 0
1386 static void cmp_disres(t_enxframe *fr1, t_enxframe *fr2, real ftol, real abstol)
1387 {
1388     int  i;
1389     char bav[64], bt[64], bs[22];
1390
1391     cmp_int(stdout, "ndisre", -1, fr1->ndisre, fr2->ndisre);
1392     if ((fr1->ndisre == fr2->ndisre) && (fr1->ndisre > 0))
1393     {
1394         sprintf(bav, "step %s: disre rav", gmx_step_str(fr1->step, bs));
1395         sprintf(bt, "step %s: disre  rt", gmx_step_str(fr1->step, bs));
1396         for (i = 0; (i < fr1->ndisre); i++)
1397         {
1398             cmp_real(stdout, bav, i, fr1->disre_rm3tav[i], fr2->disre_rm3tav[i], ftol, abstol);
1399             cmp_real(stdout, bt, i, fr1->disre_rt[i], fr2->disre_rt[i], ftol, abstol);
1400         }
1401     }
1402 }
1403 #endif
1404
1405 static void cmp_eblocks(t_enxframe *fr1, t_enxframe *fr2, real ftol, real abstol)
1406 {
1407     int  i, j, k;
1408     char buf[64], bs[22];
1409
1410     cmp_int(stdout, "nblock", -1, fr1->nblock, fr2->nblock);
1411     if ((fr1->nblock == fr2->nblock) && (fr1->nblock > 0))
1412     {
1413         for (j = 0; (j < fr1->nblock); j++)
1414         {
1415             t_enxblock *b1, *b2; /* convenience vars */
1416
1417             b1 = &(fr1->block[j]);
1418             b2 = &(fr2->block[j]);
1419
1420             sprintf(buf, "step %s: block[%d]", gmx_step_str(fr1->step, bs), j);
1421             cmp_int(stdout, buf, -1, b1->nsub, b2->nsub);
1422             cmp_int(stdout, buf, -1, b1->id, b2->id);
1423
1424             if ( (b1->nsub == b2->nsub) && (b1->id == b2->id) )
1425             {
1426                 for (i = 0; i < b1->nsub; i++)
1427                 {
1428                     t_enxsubblock *s1, *s2;
1429
1430                     s1 = &(b1->sub[i]);
1431                     s2 = &(b2->sub[i]);
1432
1433                     cmp_int(stdout, buf, -1, (int)s1->type, (int)s2->type);
1434                     cmp_int64(stdout, buf, s1->nr, s2->nr);
1435
1436                     if ((s1->type == s2->type) && (s1->nr == s2->nr))
1437                     {
1438                         switch (s1->type)
1439                         {
1440                             case xdr_datatype_float:
1441                                 for (k = 0; k < s1->nr; k++)
1442                                 {
1443                                     cmp_float(stdout, buf, i,
1444                                               s1->fval[k], s2->fval[k],
1445                                               ftol, abstol);
1446                                 }
1447                                 break;
1448                             case xdr_datatype_double:
1449                                 for (k = 0; k < s1->nr; k++)
1450                                 {
1451                                     cmp_double(stdout, buf, i,
1452                                                s1->dval[k], s2->dval[k],
1453                                                ftol, abstol);
1454                                 }
1455                                 break;
1456                             case xdr_datatype_int:
1457                                 for (k = 0; k < s1->nr; k++)
1458                                 {
1459                                     cmp_int(stdout, buf, i,
1460                                             s1->ival[k], s2->ival[k]);
1461                                 }
1462                                 break;
1463                             case xdr_datatype_int64:
1464                                 for (k = 0; k < s1->nr; k++)
1465                                 {
1466                                     cmp_int64(stdout, buf,
1467                                               s1->lval[k], s2->lval[k]);
1468                                 }
1469                                 break;
1470                             case xdr_datatype_char:
1471                                 for (k = 0; k < s1->nr; k++)
1472                                 {
1473                                     cmp_uc(stdout, buf, i,
1474                                            s1->cval[k], s2->cval[k]);
1475                                 }
1476                                 break;
1477                             case xdr_datatype_string:
1478                                 for (k = 0; k < s1->nr; k++)
1479                                 {
1480                                     cmp_str(stdout, buf, i,
1481                                             s1->sval[k], s2->sval[k]);
1482                                 }
1483                                 break;
1484                             default:
1485                                 gmx_incons("Unknown data type!!");
1486                         }
1487                     }
1488                 }
1489             }
1490         }
1491     }
1492 }
1493
1494 void comp_enx(const char *fn1, const char *fn2, real ftol, real abstol, const char *lastener)
1495 {
1496     int            nre, nre1, nre2;
1497     ener_file_t    in1, in2;
1498     int            i, j, maxener, *ind1, *ind2, *have;
1499     gmx_enxnm_t   *enm1 = nullptr, *enm2 = nullptr;
1500     t_enxframe    *fr1, *fr2;
1501     gmx_bool       b1, b2;
1502
1503     fprintf(stdout, "comparing energy file %s and %s\n\n", fn1, fn2);
1504
1505     in1 = open_enx(fn1, "r");
1506     in2 = open_enx(fn2, "r");
1507     do_enxnms(in1, &nre1, &enm1);
1508     do_enxnms(in2, &nre2, &enm2);
1509     if (nre1 != nre2)
1510     {
1511         fprintf(stdout, "There are %d and %d terms in the energy files\n\n",
1512                 nre1, nre2);
1513     }
1514     else
1515     {
1516         fprintf(stdout, "There are %d terms in the energy files\n\n", nre1);
1517     }
1518
1519     snew(ind1, nre1);
1520     snew(ind2, nre2);
1521     snew(have, nre2);
1522     nre = 0;
1523     for (i = 0; i < nre1; i++)
1524     {
1525         for (j = 0; j < nre2; j++)
1526         {
1527             if (enernm_equal(enm1[i].name, enm2[j].name))
1528             {
1529                 ind1[nre] = i;
1530                 ind2[nre] = j;
1531                 have[j]   = 1;
1532                 nre++;
1533                 break;
1534             }
1535         }
1536         if (nre == 0 || ind1[nre-1] != i)
1537         {
1538             cmp_str(stdout, "enm", i, enm1[i].name, "-");
1539         }
1540     }
1541     for (i = 0; i < nre2; i++)
1542     {
1543         if (have[i] == 0)
1544         {
1545             cmp_str(stdout, "enm", i, "-", enm2[i].name);
1546         }
1547     }
1548
1549     maxener = nre;
1550     for (i = 0; i < nre; i++)
1551     {
1552         if ((lastener != nullptr) && (std::strstr(enm1[i].name, lastener) != nullptr))
1553         {
1554             maxener = i+1;
1555             break;
1556         }
1557     }
1558
1559     fprintf(stdout, "There are %d terms to compare in the energy files\n\n",
1560             maxener);
1561
1562     for (i = 0; i < maxener; i++)
1563     {
1564         cmp_str(stdout, "unit", i, enm1[ind1[i]].unit, enm2[ind2[i]].unit);
1565     }
1566
1567     snew(fr1, 1);
1568     snew(fr2, 1);
1569     do
1570     {
1571         b1 = do_enx(in1, fr1);
1572         b2 = do_enx(in2, fr2);
1573         if (b1 && !b2)
1574         {
1575             fprintf(stdout, "\nEnd of file on %s but not on %s\n", fn2, fn1);
1576         }
1577         else if (!b1 && b2)
1578         {
1579             fprintf(stdout, "\nEnd of file on %s but not on %s\n", fn1, fn2);
1580         }
1581         else if (!b1 && !b2)
1582         {
1583             fprintf(stdout, "\nFiles read successfully\n");
1584         }
1585         else
1586         {
1587             cmp_real(stdout, "t", -1, fr1->t, fr2->t, ftol, abstol);
1588             cmp_int(stdout, "step", -1, fr1->step, fr2->step);
1589             /* We don't want to print the nre mismatch for every frame */
1590             /* cmp_int(stdout,"nre",-1,fr1->nre,fr2->nre); */
1591             if ((fr1->nre >= nre) && (fr2->nre >= nre))
1592             {
1593                 cmp_energies(stdout, fr1->step, fr1->step, fr1->ener, fr2->ener,
1594                              enm1, ftol, abstol, nre, ind1, ind2, maxener);
1595             }
1596             /*cmp_disres(fr1,fr2,ftol,abstol);*/
1597             cmp_eblocks(fr1, fr2, ftol, abstol);
1598         }
1599     }
1600     while (b1 && b2);
1601
1602     close_enx(in1);
1603     close_enx(in2);
1604
1605     free_enxframe(fr2);
1606     sfree(fr2);
1607     free_enxframe(fr1);
1608     sfree(fr1);
1609 }