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, 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, "%s\n%s", msg,
 340                   "If you want to use the correct frames before the corrupted frame and avoid this "
 341                   "fatal error set the env.var. GMX_ENX_NO_FATAL");
 342     }
 343 }
 344
 345 static void edr_strings(XDR* xdr, gmx_bool bRead, int file_version, int n, gmx_enxnm_t** nms)
 346 {
 347     int          i;
 348     gmx_enxnm_t* nm;
 349
 350     if (*nms == nullptr)
 351     {
 352         snew(*nms, n);
 353     }
 354     for (i = 0; i < n; i++)
 355     {
 356         nm = &(*nms)[i];
 357         if (bRead)
 358         {
 359             if (nm->name)
 360             {
 361                 sfree(nm->name);
 362                 nm->name = nullptr;
 363             }
 364             if (nm->unit)
 365             {
 366                 sfree(nm->unit);
 367                 nm->unit = nullptr;
 368             }
 369         }
 370         if (!xdr_string(xdr, &(nm->name), STRLEN))
 371         {
 372             gmx_file("Cannot write energy names to file; maybe you are out of disk space?");
 373         }
 374         if (file_version >= 2)
 375         {
 376             if (!xdr_string(xdr, &(nm->unit), STRLEN))
 377             {
 378                 gmx_file("Cannot write energy names to file; maybe you are out of disk space?");
 379             }
 380         }
 381         else
 382         {
 383             nm->unit = gmx_strdup("kJ/mol");
 384         }
 385     }
 386 }
 387
 388 void do_enxnms(ener_file_t ef, int* nre, gmx_enxnm_t** nms)
 389 {
 390     int      magic = -55555;
 391     XDR*     xdr;
 392     gmx_bool bRead = gmx_fio_getread(ef->fio);
 393     int      file_version;
 394
 395     xdr = gmx_fio_getxdr(ef->fio);
 396
 397     if (!xdr_int(xdr, &magic))
 398     {
 399         if (!bRead)
 400         {
 401             gmx_file("Cannot write energy names to file; maybe you are out of disk space?");
 402         }
 403         *nre = 0;
 404         return;
 405     }
 406     if (magic > 0)
 407     {
 408         /* Assume this is an old edr format */
 409         file_version          = 1;
 410         *nre                  = magic;
 411         ef->eo.bOldFileOpen   = TRUE;
 412         ef->eo.bReadFirstStep = FALSE;
 413         srenew(ef->eo.ener_prev, *nre);
 414     }
 415     else
 416     {
 417         ef->eo.bOldFileOpen = FALSE;
 418
 419         if (magic != -55555)
 420         {
 421             gmx_fatal(FARGS, "Energy names magic number mismatch, this is not a GROMACS edr file");
 422         }
 423         file_version = enx_version;
 424         xdr_int(xdr, &file_version);
 425         if (file_version > enx_version)
 426         {
 427             gmx_fatal(FARGS, "reading tpx file (%s) version %d with version %d program",
 428                       gmx_fio_getname(ef->fio), file_version, enx_version);
 429         }
 430         xdr_int(xdr, nre);
 431     }
 432     if (file_version != enx_version)
 433     {
 434         fprintf(stderr, "Note: enx file_version %d, software version %d\n", file_version, enx_version);
 435     }
 436
 437     edr_strings(xdr, bRead, file_version, *nre, nms);
 438 }
 439
 440 static gmx_bool do_eheader(ener_file_t ef,
 441                            int*        file_version,
 442                            t_enxframe* fr,
 443                            int         nre_test,
 444                            gmx_bool*   bWrongPrecision,
 445                            gmx_bool*   bOK)
 446 {
 447     int      magic = -7777777;
 448     real     first_real_to_check;
 449     int      b, zero = 0, dum = 0;
 450     gmx_bool bRead  = gmx_fio_getread(ef->fio);
 451     int      ndisre = 0;
 452     int      startb = 0;
 453 #if !GMX_DOUBLE
 454     xdr_datatype dtreal = xdr_datatype_float;
 455 #else
 456     xdr_datatype dtreal = xdr_datatype_double;
 457 #endif
 458
 459     if (bWrongPrecision)
 460     {
 461         *bWrongPrecision = FALSE;
 462     }
 463
 464     *bOK = TRUE;
 465     /* The original energy frame started with a real,
 466      * so we have to use a real for compatibility.
 467      * This is VERY DIRTY code, since do_eheader can be called
 468      * with the wrong precision set and then we could read r > -1e10,
 469      * while actually the intention was r < -1e10.
 470      * When nre_test >= 0, do_eheader should therefore terminate
 471      * before the number of i/o calls starts depending on what has been read
 472      * (which is the case for for instance the block sizes for variable
 473      * number of blocks, where this number is read before).
 474      */
 475     first_real_to_check = -2e10;
 476     if (!gmx_fio_do_real(ef->fio, first_real_to_check))
 477     {
 478         return FALSE;
 479     }
 480     if (first_real_to_check > -1e10)
 481     {
 482         /* Assume we are reading an old format */
 483         *file_version = 1;
 484         fr->t         = first_real_to_check;
 485         if (!gmx_fio_do_int(ef->fio, dum))
 486         {
 487             *bOK = FALSE;
 488         }
 489         fr->step = dum;
 490     }
 491     else
 492     {
 493         if (!gmx_fio_do_int(ef->fio, magic))
 494         {
 495             *bOK = FALSE;
 496         }
 497         if (magic != -7777777)
 498         {
 499             enx_warning("Energy header magic number mismatch, this is not a GROMACS edr file");
 500             *bOK = FALSE;
 501             return FALSE;
 502         }
 503         *file_version = enx_version;
 504         if (!gmx_fio_do_int(ef->fio, *file_version))
 505         {
 506             *bOK = FALSE;
 507         }
 508         if (*bOK && *file_version > enx_version)
 509         {
 510             gmx_fatal(FARGS, "reading tpx file (%s) version %d with version %d program",
 511                       gmx_fio_getname(ef->fio), *file_version, enx_version);
 512         }
 513         if (!gmx_fio_do_double(ef->fio, fr->t))
 514         {
 515             *bOK = FALSE;
 516         }
 517         if (!gmx_fio_do_int64(ef->fio, fr->step))
 518         {
 519             *bOK = FALSE;
 520         }
 521         if (!bRead && fr->nsum == 1)
 522         {
 523             /* Do not store sums of length 1,
 524              * since this does not add information.
 525              */
 526             if (!gmx_fio_do_int(ef->fio, zero))
 527             {
 528                 *bOK = FALSE;
 529             }
 530         }
 531         else
 532         {
 533             if (!gmx_fio_do_int(ef->fio, fr->nsum))
 534             {
 535                 *bOK = FALSE;
 536             }
 537         }
 538         if (*file_version >= 3)
 539         {
 540             if (!gmx_fio_do_int64(ef->fio, fr->nsteps))
 541             {
 542                 *bOK = FALSE;
 543             }
 544         }
 545         else
 546         {
 547             fr->nsteps = std::max(1, fr->nsum);
 548         }
 549         if (*file_version >= 5)
 550         {
 551             if (!gmx_fio_do_double(ef->fio, fr->dt))
 552             {
 553                 *bOK = FALSE;
 554             }
 555         }
 556         else
 557         {
 558             fr->dt = 0;
 559         }
 560     }
 561     if (!gmx_fio_do_int(ef->fio, fr->nre))
 562     {
 563         *bOK = FALSE;
 564     }
 565     if (*file_version < 4)
 566     {
 567         if (!gmx_fio_do_int(ef->fio, ndisre))
 568         {
 569             *bOK = FALSE;
 570         }
 571     }
 572     else
 573     {
 574         /* now reserved for possible future use */
 575         if (!gmx_fio_do_int(ef->fio, dum))
 576         {
 577             *bOK = FALSE;
 578         }
 579     }
 580
 581     if (!gmx_fio_do_int(ef->fio, fr->nblock))
 582     {
 583         *bOK = FALSE;
 584     }
 585     if (fr->nblock < 0)
 586     {
 587         *bOK = FALSE;
 588     }
 589
 590     if (ndisre != 0)
 591     {
 592         if (*file_version >= 4)
 593         {
 594             enx_warning("Distance restraint blocks in old style in new style file");
 595             *bOK = FALSE;
 596             return FALSE;
 597         }
 598         fr->nblock += 1;
 599     }
 600
 601
 602     /* Frames could have nre=0, so we can not rely only on the fr->nre check */
 603     if (bRead && nre_test >= 0
 604         && ((fr->nre > 0 && fr->nre != nre_test) || fr->nre < 0 || ndisre < 0 || fr->nblock < 0))
 605     {
 606         if (bWrongPrecision)
 607         {
 608             *bWrongPrecision = TRUE;
 609         }
 610         return *bOK;
 611     }
 612
 613     /* we now know what these should be, or we've already bailed out because
 614        of wrong precision */
 615     if (*file_version == 1 && (fr->t < 0 || fr->t > 1e20 || fr->step < 0))
 616     {
 617         enx_warning(
 618                 "edr file with negative step number or unreasonable time (and without version "
 619                 "number).");
 620         *bOK = FALSE;
 621         return FALSE;
 622     }
 623
 624
 625     if (*bOK && bRead)
 626     {
 627         add_blocks_enxframe(fr, fr->nblock);
 628     }
 629
 630     startb = 0;
 631     if (ndisre > 0)
 632     {
 633         /* sub[0] is the instantaneous data, sub[1] is time averaged */
 634         add_subblocks_enxblock(&(fr->block[0]), 2);
 635         fr->block[0].id          = enxDISRE;
 636         fr->block[0].sub[0].nr   = ndisre;
 637         fr->block[0].sub[1].nr   = ndisre;
 638         fr->block[0].sub[0].type = dtreal;
 639         fr->block[0].sub[1].type = dtreal;
 640         startb++;
 641     }
 642
 643     /* read block header info */
 644     for (b = startb; b < fr->nblock; b++)
 645     {
 646         if (*file_version < 4)
 647         {
 648             /* blocks in old version files always have 1 subblock that
 649                consists of reals. */
 650             int nrint;
 651
 652             if (bRead)
 653             {
 654                 add_subblocks_enxblock(&(fr->block[b]), 1);
 655             }
 656             else
 657             {
 658                 if (fr->block[b].nsub != 1)
 659                 {
 660                     gmx_incons("Writing an old version .edr file with too many subblocks");
 661                 }
 662                 if (fr->block[b].sub[0].type != dtreal)
 663                 {
 664                     gmx_incons("Writing an old version .edr file the wrong subblock type");
 665                 }
 666             }
 667             nrint = fr->block[b].sub[0].nr;
 668
 669             if (!gmx_fio_do_int(ef->fio, nrint))
 670             {
 671                 *bOK = FALSE;
 672             }
 673             fr->block[b].id          = b - startb;
 674             fr->block[b].sub[0].nr   = nrint;
 675             fr->block[b].sub[0].type = dtreal;
 676         }
 677         else
 678         {
 679             int i;
 680             /* in the new version files, the block header only contains
 681                the ID and the number of subblocks */
 682             int nsub = fr->block[b].nsub;
 683             *bOK     = *bOK && gmx_fio_do_int(ef->fio, fr->block[b].id);
 684             *bOK     = *bOK && gmx_fio_do_int(ef->fio, nsub);
 685
 686             fr->block[b].nsub = nsub;
 687             if (bRead)
 688             {
 689                 add_subblocks_enxblock(&(fr->block[b]), nsub);
 690             }
 691
 692             /* read/write type & size for each subblock */
 693             for (i = 0; i < nsub; i++)
 694             {
 695                 t_enxsubblock* sub    = &(fr->block[b].sub[i]); /* shortcut */
 696                 int            typenr = sub->type;
 697
 698                 *bOK = *bOK && gmx_fio_do_int(ef->fio, typenr);
 699                 *bOK = *bOK && gmx_fio_do_int(ef->fio, sub->nr);
 700
 701                 sub->type = static_cast<xdr_datatype>(typenr);
 702             }
 703         }
 704     }
 705     if (!gmx_fio_do_int(ef->fio, fr->e_size))
 706     {
 707         *bOK = FALSE;
 708     }
 709
 710     /* now reserved for possible future use */
 711     if (!gmx_fio_do_int(ef->fio, dum))
 712     {
 713         *bOK = FALSE;
 714     }
 715
 716     /* Do a dummy int to keep the format compatible with the old code */
 717     if (!gmx_fio_do_int(ef->fio, dum))
 718     {
 719         *bOK = FALSE;
 720     }
 721
 722     if (*bOK && *file_version == 1 && nre_test < 0)
 723     {
 724         if (!ef->eo.bReadFirstStep)
 725         {
 726             ef->eo.bReadFirstStep = TRUE;
 727             ef->eo.first_step     = fr->step;
 728             ef->eo.step_prev      = fr->step;
 729             ef->eo.nsum_prev      = 0;
 730         }
 731
 732         fr->nsum   = fr->step - ef->eo.first_step + 1;
 733         fr->nsteps = fr->step - ef->eo.step_prev;
 734         fr->dt     = 0;
 735     }
 736
 737     return *bOK;
 738 }
 739
 740 void free_enxnms(int n, gmx_enxnm_t* nms)
 741 {
 742     int i;
 743
 744     for (i = 0; i < n; i++)
 745     {
 746         sfree(nms[i].name);
 747         sfree(nms[i].unit);
 748     }
 749
 750     sfree(nms);
 751 }
 752
 753 void close_enx(ener_file_t ef)
 754 {
 755     if (ef == nullptr)
 756     {
 757         // Nothing to do
 758         return;
 759     }
 760     if (gmx_fio_close(ef->fio) != 0)
 761     {
 762         gmx_file(
 763                 "Cannot close energy file; it might be corrupt, or maybe you are out of disk "
 764                 "space?");
 765     }
 766 }
 767
 768 void done_ener_file(ener_file_t ef)
 769 {
 770     // Free the contents, then the pointer itself
 771     close_enx(ef);
 772     sfree(ef);
 773 }
 774
 775 /*!\brief Return TRUE if a file exists but is empty, otherwise FALSE.
 776  *
 777  * If the file exists but has length larger than zero, if it does not exist, or
 778  * if there is a file system error, FALSE will be returned instead.
 779  *
 780  * \param fn  File name to test
 781  *
 782  * \return TRUE if file could be open but is empty, otherwise FALSE.
 783  */
 784 static gmx_bool empty_file(const char* fn)
 785 {
 786     FILE*    fp;
 787     char     dum;
 788     int      ret;
 789     gmx_bool bEmpty;
 790
 791     fp     = gmx_fio_fopen(fn, "r");
 792     ret    = fread(&dum, sizeof(dum), 1, fp);
 793     bEmpty = (feof(fp) != 0);
 794     gmx_fio_fclose(fp);
 795
 796     // bEmpty==TRUE but ret!=0 would likely be some strange I/O error, but at
 797     // least it's not a normal empty file, so we return FALSE in that case.
 798     return (bEmpty && ret == 0);
 799 }
 800
 801
 802 ener_file_t open_enx(const char* fn, const char* mode)
 803 {
 804     int               nre;
 805     gmx_enxnm_t*      nms          = nullptr;
 806     int               file_version = -1;
 807     t_enxframe*       fr;
 808     gmx_bool          bWrongPrecision, bOK = TRUE;
 809     struct ener_file* ef;
 810
 811     snew(ef, 1);
 812
 813     if (mode[0] == 'r')
 814     {
 815         ef->fio = gmx_fio_open(fn, mode);
 816         gmx_fio_setprecision(ef->fio, FALSE);
 817         do_enxnms(ef, &nre, &nms);
 818         snew(fr, 1);
 819         do_eheader(ef, &file_version, fr, nre, &bWrongPrecision, &bOK);
 820         if (!bOK)
 821         {
 822             gmx_file("Cannot read energy file header. Corrupt file?");
 823         }
 824
 825         /* Now check whether this file is in single precision */
 826         if (!bWrongPrecision
 827             && ((fr->e_size && (fr->nre == nre)
 828                  && (nre * 4 * static_cast<long int>(sizeof(float)) == fr->e_size))))
 829         {
 830             fprintf(stderr, "Opened %s as single precision energy file\n", fn);
 831             free_enxnms(nre, nms);
 832         }
 833         else
 834         {
 835             gmx_fio_rewind(ef->fio);
 836             gmx_fio_setprecision(ef->fio, TRUE);
 837             do_enxnms(ef, &nre, &nms);
 838             do_eheader(ef, &file_version, fr, nre, &bWrongPrecision, &bOK);
 839             if (!bOK)
 840             {
 841                 gmx_file("Cannot write energy file header; maybe you are out of disk space?");
 842             }
 843
 844             if (((fr->e_size && (fr->nre == nre)
 845                   && (nre * 4 * static_cast<long int>(sizeof(double)) == fr->e_size))))
 846             {
 847                 fprintf(stderr, "Opened %s as double precision energy file\n", fn);
 848             }
 849             else
 850             {
 851                 if (empty_file(fn))
 852                 {
 853                     gmx_fatal(FARGS, "File %s is empty", fn);
 854                 }
 855                 else
 856                 {
 857                     gmx_fatal(FARGS, "Energy file %s not recognized, maybe different CPU?", fn);
 858                 }
 859             }
 860             free_enxnms(nre, nms);
 861         }
 862         free_enxframe(fr);
 863         sfree(fr);
 864         gmx_fio_rewind(ef->fio);
 865     }
 866     else
 867     {
 868         ef->fio = gmx_fio_open(fn, mode);
 869     }
 870
 871     ef->framenr   = 0;
 872     ef->frametime = 0;
 873     return ef;
 874 }
 875
 876 t_fileio* enx_file_pointer(const ener_file* ef)
 877 {
 878     return ef->fio;
 879 }
 880
 881 static void convert_full_sums(ener_old_t* ener_old, t_enxframe* fr)
 882 {
 883     int    nstep_all;
 884     int    ne, ns, i;
 885     double esum_all, eav_all;
 886
 887     if (fr->nsum > 0)
 888     {
 889         ne = 0;
 890         ns = 0;
 891         for (i = 0; i < fr->nre; i++)
 892         {
 893             if (fr->ener[i].e != 0)
 894             {
 895                 ne++;
 896             }
 897             if (fr->ener[i].esum != 0)
 898             {
 899                 ns++;
 900             }
 901         }
 902         if (ne > 0 && ns == 0)
 903         {
 904             /* We do not have all energy sums */
 905             fr->nsum = 0;
 906         }
 907     }
 908
 909     /* Convert old full simulation sums to sums between energy frames */
 910     nstep_all = fr->step - ener_old->first_step + 1;
 911     if (fr->nsum > 1 && fr->nsum == nstep_all && ener_old->nsum_prev > 0)
 912     {
 913         /* Set the new sum length: the frame step difference */
 914         fr->nsum = fr->step - ener_old->step_prev;
 915         for (i = 0; i < fr->nre; i++)
 916         {
 917             esum_all         = fr->ener[i].esum;
 918             eav_all          = fr->ener[i].eav;
 919             fr->ener[i].esum = esum_all - ener_old->ener_prev[i].esum;
 920             fr->ener[i].eav =
 921                     eav_all - ener_old->ener_prev[i].eav
 922                     - gmx::square(ener_old->ener_prev[i].esum / (nstep_all - fr->nsum) - esum_all / nstep_all)
 923                               * (nstep_all - fr->nsum) * nstep_all / static_cast<double>(fr->nsum);
 924             ener_old->ener_prev[i].esum = esum_all;
 925             ener_old->ener_prev[i].eav  = eav_all;
 926         }
 927         ener_old->nsum_prev = nstep_all;
 928     }
 929     else if (fr->nsum > 0)
 930     {
 931         if (fr->nsum != nstep_all)
 932         {
 933             fprintf(stderr,
 934                     "\nWARNING: something is wrong with the energy sums, will not use exact "
 935                     "averages\n");
 936             ener_old->nsum_prev = 0;
 937         }
 938         else
 939         {
 940             ener_old->nsum_prev = nstep_all;
 941         }
 942         /* Copy all sums to ener_prev */
 943         for (i = 0; i < fr->nre; i++)
 944         {
 945             ener_old->ener_prev[i].esum = fr->ener[i].esum;
 946             ener_old->ener_prev[i].eav  = fr->ener[i].eav;
 947         }
 948     }
 949
 950     ener_old->step_prev = fr->step;
 951 }
 952
 953 gmx_bool do_enx(ener_file_t ef, t_enxframe* fr)
 954 {
 955     int      file_version = -1;
 956     int      i, b;
 957     gmx_bool bRead, bOK, bOK1, bSane;
 958     real     tmp1, tmp2, rdum;
 959     /*int       d_size;*/
 960
 961     bOK   = TRUE;
 962     bRead = gmx_fio_getread(ef->fio);
 963     if (!bRead)
 964     {
 965         fr->e_size = fr->nre * sizeof(fr->ener[0].e) * 4;
 966         /*d_size = fr->ndisre*(sizeof(real)*2);*/
 967     }
 968
 969     if (!do_eheader(ef, &file_version, fr, -1, nullptr, &bOK))
 970     {
 971         if (bRead)
 972         {
 973             fprintf(stderr, "\rLast energy frame read %d time %8.3f         ", ef->framenr - 1,
 974                     ef->frametime);
 975             fflush(stderr);
 976
 977             if (!bOK)
 978             {
 979                 fprintf(stderr, "\nWARNING: Incomplete energy frame: nr %d time %8.3f\n",
 980                         ef->framenr, fr->t);
 981             }
 982         }
 983         else
 984         {
 985             gmx_file("Cannot write energy file header; maybe you are out of disk space?");
 986         }
 987         return FALSE;
 988     }
 989     if (bRead)
 990     {
 991         if ((ef->framenr < 20 || ef->framenr % 10 == 0) && (ef->framenr < 200 || ef->framenr % 100 == 0)
 992             && (ef->framenr < 2000 || ef->framenr % 1000 == 0))
 993         {
 994             fprintf(stderr, "\rReading energy frame %6d time %8.3f         ", ef->framenr, fr->t);
 995         }
 996         ef->framenr++;
 997         ef->frametime = fr->t;
 998     }
 999     /* Check sanity of this header */
1000     bSane = fr->nre > 0;
1001     for (b = 0; b < fr->nblock; b++)
1002     {
1003         bSane = bSane || (fr->block[b].nsub > 0);
1004     }
1005     if (!((fr->step >= 0) && bSane) && bRead)
1006     {
1007         fprintf(stderr, "\nWARNING: there may be something wrong with energy file %s\n",
1008                 gmx_fio_getname(ef->fio));
1009         fprintf(stderr, "Found: step=%" PRId64 ", nre=%d, nblock=%d, time=%g.\n", fr->step, fr->nre,
1010                 fr->nblock, fr->t);
1011     }
1012     if (bRead && fr->nre > fr->e_alloc)
1013     {
1014         srenew(fr->ener, fr->nre);
1015         for (i = fr->e_alloc; (i < fr->nre); i++)
1016         {
1017             fr->ener[i].e    = 0;
1018             fr->ener[i].eav  = 0;
1019             fr->ener[i].esum = 0;
1020         }
1021         fr->e_alloc = fr->nre;
1022     }
1023
1024     for (i = 0; i < fr->nre; i++)
1025     {
1026         bOK = bOK && gmx_fio_do_real(ef->fio, fr->ener[i].e);
1027
1028         /* Do not store sums of length 1,
1029          * since this does not add information.
1030          */
1031         if (file_version == 1 || (bRead && fr->nsum > 0) || fr->nsum > 1)
1032         {
1033             tmp1 = fr->ener[i].eav;
1034             bOK  = bOK && gmx_fio_do_real(ef->fio, tmp1);
1035             if (bRead)
1036             {
1037                 fr->ener[i].eav = tmp1;
1038             }
1039
1040             /* This is to save only in single precision (unless compiled in DP) */
1041             tmp2 = fr->ener[i].esum;
1042             bOK  = bOK && gmx_fio_do_real(ef->fio, tmp2);
1043             if (bRead)
1044             {
1045                 fr->ener[i].esum = tmp2;
1046             }
1047
1048             if (file_version == 1)
1049             {
1050                 /* Old, unused real */
1051                 rdum = 0;
1052                 bOK  = bOK && gmx_fio_do_real(ef->fio, rdum);
1053             }
1054         }
1055     }
1056
1057     /* Here we can not check for file_version==1, since one could have
1058      * continued an old format simulation with a new one with mdrun -append.
1059      */
1060     if (bRead && ef->eo.bOldFileOpen)
1061     {
1062         /* Convert old full simulation sums to sums between energy frames */
1063         convert_full_sums(&(ef->eo), fr);
1064     }
1065     /* read the blocks */
1066     for (b = 0; b < fr->nblock; b++)
1067     {
1068         /* now read the subblocks. */
1069         int nsub = fr->block[b].nsub; /* shortcut */
1070         int i;
1071
1072         for (i = 0; i < nsub; i++)
1073         {
1074             t_enxsubblock* sub = &(fr->block[b].sub[i]); /* shortcut */
1075
1076             if (bRead)
1077             {
1078                 enxsubblock_alloc(sub);
1079             }
1080
1081             /* read/write data */
1082             switch (sub->type)
1083             {
1084                 case xdr_datatype_float:
1085                     bOK1 = gmx_fio_ndo_float(ef->fio, sub->fval, sub->nr);
1086                     break;
1087                 case xdr_datatype_double:
1088                     bOK1 = gmx_fio_ndo_double(ef->fio, sub->dval, sub->nr);
1089                     break;
1090                 case xdr_datatype_int: bOK1 = gmx_fio_ndo_int(ef->fio, sub->ival, sub->nr); break;
1091                 case xdr_datatype_int64:
1092                     bOK1 = gmx_fio_ndo_int64(ef->fio, sub->lval, sub->nr);
1093                     break;
1094                 case xdr_datatype_char:
1095                     bOK1 = gmx_fio_ndo_uchar(ef->fio, sub->cval, sub->nr);
1096                     break;
1097                 case xdr_datatype_string:
1098                     bOK1 = gmx_fio_ndo_string(ef->fio, sub->sval, sub->nr);
1099                     break;
1100                 default:
1101                     gmx_incons(
1102                             "Reading unknown block data type: this file is corrupted or from the "
1103                             "future");
1104             }
1105             bOK = bOK && bOK1;
1106         }
1107     }
1108
1109     if (!bRead)
1110     {
1111         if (gmx_fio_flush(ef->fio) != 0)
1112         {
1113             gmx_file("Cannot write energy file; maybe you are out of disk space?");
1114         }
1115     }
1116
1117     if (!bOK)
1118     {
1119         if (bRead)
1120         {
1121             fprintf(stderr, "\nLast energy frame read %d", ef->framenr - 1);
1122             fprintf(stderr, "\nWARNING: Incomplete energy frame: nr %d time %8.3f\n", ef->framenr, fr->t);
1123         }
1124         else
1125         {
1126             gmx_fatal(FARGS, "could not write energies");
1127         }
1128         return FALSE;
1129     }
1130
1131     return TRUE;
1132 }
1133
1134 static real find_energy(const char* name, int nre, gmx_enxnm_t* enm, t_enxframe* fr)
1135 {
1136     int i;
1137
1138     for (i = 0; i < nre; i++)
1139     {
1140         if (std::strcmp(enm[i].name, name) == 0)
1141         {
1142             return fr->ener[i].e;
1143         }
1144     }
1145
1146     gmx_fatal(FARGS,
1147               "Could not find energy term named '%s'. Either the energy file is from a different "
1148               "run or this state variable is not stored in the energy file. In the latter case "
1149               "(and if you did not modify the T/P-coupling setup), you can read the state in mdrun "
1150               "instead, by passing in a checkpoint file.",
1151               name);
1152 }
1153
1154
1155 void get_enx_state(const char* fn, real t, const SimulationGroups& groups, t_inputrec* ir, t_state* state)
1156 {
1157     /* Should match the names in mdebin.c */
1158     static const char* boxvel_nm[] = { "Box-Vel-XX", "Box-Vel-YY", "Box-Vel-ZZ",
1159                                        "Box-Vel-YX", "Box-Vel-ZX", "Box-Vel-ZY" };
1160
1161     static const char* baro_nm[] = { "Barostat" };
1162
1163
1164     int          ind0[] = { XX, YY, ZZ, YY, ZZ, ZZ };
1165     int          ind1[] = { XX, YY, ZZ, XX, XX, YY };
1166     int          nre, nfr, i, j, ni, npcoupl;
1167     char         buf[STRLEN];
1168     const char*  bufi;
1169     gmx_enxnm_t* enm = nullptr;
1170     t_enxframe*  fr;
1171     ener_file_t  in;
1172
1173     in = open_enx(fn, "r");
1174     do_enxnms(in, &nre, &enm);
1175     snew(fr, 1);
1176     nfr = 0;
1177     while ((nfr == 0 || fr->t != t) && do_enx(in, fr))
1178     {
1179         nfr++;
1180     }
1181     close_enx(in);
1182     fprintf(stderr, "\n");
1183
1184     if (nfr == 0 || fr->t != t)
1185     {
1186         gmx_fatal(FARGS, "Could not find frame with time %f in '%s'", t, fn);
1187     }
1188
1189     npcoupl = TRICLINIC(ir->compress) ? 6 : 3;
1190     if (ir->epc == epcPARRINELLORAHMAN)
1191     {
1192         clear_mat(state->boxv);
1193         for (i = 0; i < npcoupl; i++)
1194         {
1195             state->boxv[ind0[i]][ind1[i]] = find_energy(boxvel_nm[i], nre, enm, fr);
1196         }
1197         fprintf(stderr, "\nREAD %d BOX VELOCITIES FROM %s\n\n", npcoupl, fn);
1198     }
1199
1200     if (ir->etc == etcNOSEHOOVER)
1201     {
1202         char cns[20];
1203
1204         cns[0] = '\0';
1205
1206         for (i = 0; i < state->ngtc; i++)
1207         {
1208             ni   = groups.groups[SimulationAtomGroupType::TemperatureCoupling][i];
1209             bufi = *(groups.groupNames[ni]);
1210             for (j = 0; (j < state->nhchainlength); j++)
1211             {
1212                 if (inputrecNvtTrotter(ir))
1213                 {
1214                     sprintf(cns, "-%d", j);
1215                 }
1216                 sprintf(buf, "Xi%s-%s", cns, bufi);
1217                 state->nosehoover_xi[i] = find_energy(buf, nre, enm, fr);
1218                 sprintf(buf, "vXi%s-%s", cns, bufi);
1219                 state->nosehoover_vxi[i] = find_energy(buf, nre, enm, fr);
1220             }
1221         }
1222         fprintf(stderr, "\nREAD %d NOSE-HOOVER Xi chains FROM %s\n\n", state->ngtc, fn);
1223
1224         if (inputrecNptTrotter(ir) || inputrecNphTrotter(ir))
1225         {
1226             for (i = 0; i < state->nnhpres; i++)
1227             {
1228                 bufi = baro_nm[0]; /* All barostat DOF's together for now */
1229                 for (j = 0; (j < state->nhchainlength); j++)
1230                 {
1231                     sprintf(buf, "Xi-%d-%s", j, bufi);
1232                     state->nhpres_xi[i] = find_energy(buf, nre, enm, fr);
1233                     sprintf(buf, "vXi-%d-%s", j, bufi);
1234                     state->nhpres_vxi[i] = find_energy(buf, nre, enm, fr);
1235                 }
1236             }
1237             fprintf(stderr, "\nREAD %d NOSE-HOOVER BAROSTAT Xi chains FROM %s\n\n", state->nnhpres, fn);
1238         }
1239     }
1240
1241     free_enxnms(nre, enm);
1242     free_enxframe(fr);
1243     sfree(fr);
1244 }
1245
1246 static real ener_tensor_diag(int          n,
1247                              const int*   ind1,
1248                              const int*   ind2,
1249                              gmx_enxnm_t* enm1,
1250                              const int*   tensi,
1251                              int          i,
1252                              t_energy     e1[],
1253                              t_energy     e2[])
1254 {
1255     int    d1, d2;
1256     int    j;
1257     real   prod1, prod2;
1258     int    nfound;
1259     size_t len;
1260
1261     d1 = tensi[i] / DIM;
1262     d2 = tensi[i] - d1 * DIM;
1263
1264     /* Find the diagonal elements d1 and d2 */
1265     len    = std::strlen(enm1[ind1[i]].name);
1266     prod1  = 1;
1267     prod2  = 1;
1268     nfound = 0;
1269     for (j = 0; j < n; j++)
1270     {
1271         if (tensi[j] >= 0 && std::strlen(enm1[ind1[j]].name) == len
1272             && std::strncmp(enm1[ind1[i]].name, enm1[ind1[j]].name, len - 2) == 0
1273             && (tensi[j] == d1 * DIM + d1 || tensi[j] == d2 * DIM + d2))
1274         {
1275             prod1 *= fabs(e1[ind1[j]].e);
1276             prod2 *= fabs(e2[ind2[j]].e);
1277             nfound++;
1278         }
1279     }
1280
1281     if (nfound == 2)
1282     {
1283         return 0.5 * (std::sqrt(prod1) + std::sqrt(prod2));
1284     }
1285     else
1286     {
1287         return 0;
1288     }
1289 }
1290
1291 static gmx_bool enernm_equal(const char* nm1, const char* nm2)
1292 {
1293     int len1, len2;
1294
1295     len1 = std::strlen(nm1);
1296     len2 = std::strlen(nm2);
1297
1298     /* Remove " (bar)" at the end of a name */
1299     if (len1 > 6 && std::strcmp(nm1 + len1 - 6, " (bar)") == 0)
1300     {
1301         len1 -= 6;
1302     }
1303     if (len2 > 6 && std::strcmp(nm2 + len2 - 6, " (bar)") == 0)
1304     {
1305         len2 -= 6;
1306     }
1307
1308     return (len1 == len2 && gmx_strncasecmp(nm1, nm2, len1) == 0);
1309 }
1310
1311 static void cmp_energies(FILE*        fp,
1312                          int          step1,
1313                          int          step2,
1314                          t_energy     e1[],
1315                          t_energy     e2[],
1316                          gmx_enxnm_t* enm1,
1317                          real         ftol,
1318                          real         abstol,
1319                          int          nre,
1320                          int*         ind1,
1321                          int*         ind2,
1322                          int          maxener)
1323 {
1324     int  i, ii;
1325     int *tensi, len, d1, d2;
1326     real ftol_i, abstol_i;
1327
1328     snew(tensi, maxener);
1329     /* Check for tensor elements ending on "-XX", "-XY", ... , "-ZZ" */
1330     for (i = 0; (i < maxener); i++)
1331     {
1332         ii       = ind1[i];
1333         tensi[i] = -1;
1334         len      = std::strlen(enm1[ii].name);
1335         if (len > 3 && enm1[ii].name[len - 3] == '-')
1336         {
1337             d1 = enm1[ii].name[len - 2] - 'X';
1338             d2 = enm1[ii].name[len - 1] - 'X';
1339             if (d1 >= 0 && d1 < DIM && d2 >= 0 && d2 < DIM)
1340             {
1341                 tensi[i] = d1 * DIM + d2;
1342             }
1343         }
1344     }
1345
1346     for (i = 0; (i < maxener); i++)
1347     {
1348         /* Check if this is an off-diagonal tensor element */
1349         if (tensi[i] >= 0 && tensi[i] != 0 && tensi[i] != 4 && tensi[i] != 8)
1350         {
1351             /* Turn on the relative tolerance check (4 is maximum relative diff.) */
1352             ftol_i = 5;
1353             /* Do the relative tolerance through an absolute tolerance times
1354              * the size of diagonal components of the tensor.
1355              */
1356             abstol_i = ftol * ener_tensor_diag(nre, ind1, ind2, enm1, tensi, i, e1, e2);
1357             if (debug)
1358             {
1359                 fprintf(debug, "tensor '%s' val %f diag %f\n", enm1[i].name, e1[i].e, abstol_i / ftol);
1360             }
1361             if (abstol_i > 0)
1362             {
1363                 /* We found a diagonal, we need to check with the minimum tolerance */
1364                 abstol_i = std::min(abstol_i, abstol);
1365             }
1366             else
1367             {
1368                 /* We did not find a diagonal, ignore the relative tolerance check */
1369                 abstol_i = abstol;
1370             }
1371         }
1372         else
1373         {
1374             ftol_i   = ftol;
1375             abstol_i = abstol;
1376         }
1377         if (!equal_real(e1[ind1[i]].e, e2[ind2[i]].e, ftol_i, abstol_i))
1378         {
1379             fprintf(fp, "%-15s  step %3d:  %12g,  step %3d: %12g\n", enm1[ind1[i]].name, step1,
1380                     e1[ind1[i]].e, step2, e2[ind2[i]].e);
1381         }
1382     }
1383
1384     sfree(tensi);
1385 }
1386
1387 #if 0
1388 static void cmp_disres(t_enxframe *fr1, t_enxframe *fr2, real ftol, real abstol)
1389 {
1390     int  i;
1391     char bav[64], bt[64], bs[22];
1392
1393     cmp_int(stdout, "ndisre", -1, fr1->ndisre, fr2->ndisre);
1394     if ((fr1->ndisre == fr2->ndisre) && (fr1->ndisre > 0))
1395     {
1396         sprintf(bav, "step %s: disre rav", gmx_step_str(fr1->step, bs));
1397         sprintf(bt, "step %s: disre  rt", gmx_step_str(fr1->step, bs));
1398         for (i = 0; (i < fr1->ndisre); i++)
1399         {
1400             cmp_real(stdout, bav, i, fr1->disre_rm3tav[i], fr2->disre_rm3tav[i], ftol, abstol);
1401             cmp_real(stdout, bt, i, fr1->disre_rt[i], fr2->disre_rt[i], ftol, abstol);
1402         }
1403     }
1404 }
1405 #endif
1406
1407 static void cmp_eblocks(t_enxframe* fr1, t_enxframe* fr2, real ftol, real abstol)
1408 {
1409     int  i, j, k;
1410     char buf[64], bs[22];
1411
1412     cmp_int(stdout, "nblock", -1, fr1->nblock, fr2->nblock);
1413     if ((fr1->nblock == fr2->nblock) && (fr1->nblock > 0))
1414     {
1415         for (j = 0; (j < fr1->nblock); j++)
1416         {
1417             t_enxblock *b1, *b2; /* convenience vars */
1418
1419             b1 = &(fr1->block[j]);
1420             b2 = &(fr2->block[j]);
1421
1422             sprintf(buf, "step %s: block[%d]", gmx_step_str(fr1->step, bs), j);
1423             cmp_int(stdout, buf, -1, b1->nsub, b2->nsub);
1424             cmp_int(stdout, buf, -1, b1->id, b2->id);
1425
1426             if ((b1->nsub == b2->nsub) && (b1->id == b2->id))
1427             {
1428                 for (i = 0; i < b1->nsub; i++)
1429                 {
1430                     t_enxsubblock *s1, *s2;
1431
1432                     s1 = &(b1->sub[i]);
1433                     s2 = &(b2->sub[i]);
1434
1435                     cmp_int(stdout, buf, -1, static_cast<int>(s1->type), static_cast<int>(s2->type));
1436                     cmp_int64(stdout, buf, s1->nr, s2->nr);
1437
1438                     if ((s1->type == s2->type) && (s1->nr == s2->nr))
1439                     {
1440                         switch (s1->type)
1441                         {
1442                             case xdr_datatype_float:
1443                                 for (k = 0; k < s1->nr; k++)
1444                                 {
1445                                     cmp_float(stdout, buf, i, s1->fval[k], s2->fval[k], ftol, abstol);
1446                                 }
1447                                 break;
1448                             case xdr_datatype_double:
1449                                 for (k = 0; k < s1->nr; k++)
1450                                 {
1451                                     cmp_double(stdout, buf, i, s1->dval[k], s2->dval[k], ftol, abstol);
1452                                 }
1453                                 break;
1454                             case xdr_datatype_int:
1455                                 for (k = 0; k < s1->nr; k++)
1456                                 {
1457                                     cmp_int(stdout, buf, i, s1->ival[k], s2->ival[k]);
1458                                 }
1459                                 break;
1460                             case xdr_datatype_int64:
1461                                 for (k = 0; k < s1->nr; k++)
1462                                 {
1463                                     cmp_int64(stdout, buf, s1->lval[k], s2->lval[k]);
1464                                 }
1465                                 break;
1466                             case xdr_datatype_char:
1467                                 for (k = 0; k < s1->nr; k++)
1468                                 {
1469                                     cmp_uc(stdout, buf, i, s1->cval[k], s2->cval[k]);
1470                                 }
1471                                 break;
1472                             case xdr_datatype_string:
1473                                 for (k = 0; k < s1->nr; k++)
1474                                 {
1475                                     cmp_str(stdout, buf, i, s1->sval[k], s2->sval[k]);
1476                                 }
1477                                 break;
1478                             default: gmx_incons("Unknown data type!!");
1479                         }
1480                     }
1481                 }
1482             }
1483         }
1484     }
1485 }
1486
1487 void comp_enx(const char* fn1, const char* fn2, real ftol, real abstol, const char* lastener)
1488 {
1489     int          nre, nre1, nre2;
1490     ener_file_t  in1, in2;
1491     int          i, j, maxener, *ind1, *ind2, *have;
1492     gmx_enxnm_t *enm1 = nullptr, *enm2 = nullptr;
1493     t_enxframe * fr1, *fr2;
1494     gmx_bool     b1, b2;
1495
1496     fprintf(stdout, "comparing energy file %s and %s\n\n", fn1, fn2);
1497
1498     in1 = open_enx(fn1, "r");
1499     in2 = open_enx(fn2, "r");
1500     do_enxnms(in1, &nre1, &enm1);
1501     do_enxnms(in2, &nre2, &enm2);
1502     if (nre1 != nre2)
1503     {
1504         fprintf(stdout, "There are %d and %d terms in the energy files\n\n", nre1, nre2);
1505     }
1506     else
1507     {
1508         fprintf(stdout, "There are %d terms in the energy files\n\n", nre1);
1509     }
1510
1511     snew(ind1, nre1);
1512     snew(ind2, nre2);
1513     snew(have, nre2);
1514     nre = 0;
1515     for (i = 0; i < nre1; i++)
1516     {
1517         for (j = 0; j < nre2; j++)
1518         {
1519             if (enernm_equal(enm1[i].name, enm2[j].name))
1520             {
1521                 ind1[nre] = i;
1522                 ind2[nre] = j;
1523                 have[j]   = 1;
1524                 nre++;
1525                 break;
1526             }
1527         }
1528         if (nre == 0 || ind1[nre - 1] != i)
1529         {
1530             cmp_str(stdout, "enm", i, enm1[i].name, "-");
1531         }
1532     }
1533     for (i = 0; i < nre2; i++)
1534     {
1535         if (have[i] == 0)
1536         {
1537             cmp_str(stdout, "enm", i, "-", enm2[i].name);
1538         }
1539     }
1540
1541     maxener = nre;
1542     for (i = 0; i < nre; i++)
1543     {
1544         if ((lastener != nullptr) && (std::strstr(enm1[i].name, lastener) != nullptr))
1545         {
1546             maxener = i + 1;
1547             break;
1548         }
1549     }
1550
1551     fprintf(stdout, "There are %d terms to compare in the energy files\n\n", maxener);
1552
1553     for (i = 0; i < maxener; i++)
1554     {
1555         cmp_str(stdout, "unit", i, enm1[ind1[i]].unit, enm2[ind2[i]].unit);
1556     }
1557
1558     snew(fr1, 1);
1559     snew(fr2, 1);
1560     do
1561     {
1562         b1 = do_enx(in1, fr1);
1563         b2 = do_enx(in2, fr2);
1564         if (b1 && !b2)
1565         {
1566             fprintf(stdout, "\nEnd of file on %s but not on %s\n", fn2, fn1);
1567         }
1568         else if (!b1 && b2)
1569         {
1570             fprintf(stdout, "\nEnd of file on %s but not on %s\n", fn1, fn2);
1571         }
1572         else if (!b1 && !b2)
1573         {
1574             fprintf(stdout, "\nFiles read successfully\n");
1575         }
1576         else
1577         {
1578             cmp_real(stdout, "t", -1, fr1->t, fr2->t, ftol, abstol);
1579             cmp_int(stdout, "step", -1, fr1->step, fr2->step);
1580             /* We don't want to print the nre mismatch for every frame */
1581             /* cmp_int(stdout,"nre",-1,fr1->nre,fr2->nre); */
1582             if ((fr1->nre >= nre) && (fr2->nre >= nre))
1583             {
1584                 cmp_energies(stdout, fr1->step, fr1->step, fr1->ener, fr2->ener, enm1, ftol, abstol,
1585                              nre, ind1, ind2, maxener);
1586             }
1587             /*cmp_disres(fr1,fr2,ftol,abstol);*/
1588             cmp_eblocks(fr1, fr2, ftol, abstol);
1589         }
1590     } while (b1 && b2);
1591
1592     close_enx(in1);
1593     close_enx(in2);
1594
1595     free_enxframe(fr2);
1596     sfree(fr2);
1597     free_enxframe(fr1);
1598     sfree(fr1);
1599 }