src/gromacs/fileio/enxio.cpp

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