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