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