src/gromacs/fileio/enxio.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2016,2017, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include "enxio.h"
  40
  41 #include <cstdlib>
  42 #include <cstring>
  43
  44 #include <algorithm>
  45
  46 #include "gromacs/fileio/gmxfio.h"
  47 #include "gromacs/fileio/gmxfio-xdr.h"
  48 #include "gromacs/fileio/xdrf.h"
  49 #include "gromacs/math/functions.h"
  50 #include "gromacs/math/vec.h"
  51 #include "gromacs/mdtypes/inputrec.h"
  52 #include "gromacs/mdtypes/md_enums.h"
  53 #include "gromacs/mdtypes/state.h"
  54 #include "gromacs/pbcutil/pbc.h"
  55 #include "gromacs/topology/topology.h"
  56 #include "gromacs/utility/compare.h"
  57 #include "gromacs/utility/fatalerror.h"
  58 #include "gromacs/utility/futil.h"
  59 #include "gromacs/utility/gmxassert.h"
  60 #include "gromacs/utility/smalloc.h"
  61
  62 /* The source code in this file should be thread-safe.
  63          Please keep it that way. */
  64
  65 /* This number should be increased whenever the file format changes! */
  66 static const int enx_version = 5;
  67
  68 const char      *enx_block_id_name[] = {
  69     "Averaged orientation restraints",
  70     "Instantaneous orientation restraints",
  71     "Orientation restraint order tensor(s)",
  72     "Distance restraints",
  73     "Free energy data",
  74     "BAR histogram",
  75     "Delta H raw data",
  76     "AWH data"
  77 };
  78
  79
  80 /* Stuff for reading pre 4.1 energy files */
  81 typedef struct {
  82     gmx_bool     bOldFileOpen;   /* Is this an open old file? */
  83     gmx_bool     bReadFirstStep; /* Did we read the first step? */
  84     int          first_step;     /* First step in the energy file */
  85     int          step_prev;      /* Previous step */
  86     int          nsum_prev;      /* Previous step sum length */
  87     t_energy    *ener_prev;      /* Previous energy sums */
  88 } ener_old_t;
  89
  90 struct ener_file
  91 {
  92     ener_old_t eo;
  93     t_fileio  *fio;
  94     int        framenr;
  95     real       frametime;
  96 };
  97
  98 static void enxsubblock_init(t_enxsubblock *sb)
  99 {
 100     sb->nr = 0;
 101 #if GMX_DOUBLE
 102     sb->type = xdr_datatype_double;
 103 #else
 104     sb->type = xdr_datatype_float;
 105 #endif
 106     sb->fval       = nullptr;
 107     sb->dval       = nullptr;
 108     sb->ival       = nullptr;
 109     sb->lval       = nullptr;
 110     sb->cval       = nullptr;
 111     sb->sval       = nullptr;
 112     sb->fval_alloc = 0;
 113     sb->dval_alloc = 0;
 114     sb->ival_alloc = 0;
 115     sb->lval_alloc = 0;
 116     sb->cval_alloc = 0;
 117     sb->sval_alloc = 0;
 118 }
 119
 120 static void enxsubblock_free(t_enxsubblock *sb)
 121 {
 122     if (sb->fval_alloc)
 123     {
 124         sfree(sb->fval);
 125         sb->fval_alloc = 0;
 126         sb->fval       = nullptr;
 127     }
 128     if (sb->dval_alloc)
 129     {
 130         sfree(sb->dval);
 131         sb->dval_alloc = 0;
 132         sb->dval       = nullptr;
 133     }
 134     if (sb->ival_alloc)
 135     {
 136         sfree(sb->ival);
 137         sb->ival_alloc = 0;
 138         sb->ival       = nullptr;
 139     }
 140     if (sb->lval_alloc)
 141     {
 142         sfree(sb->lval);
 143         sb->lval_alloc = 0;
 144         sb->lval       = nullptr;
 145     }
 146     if (sb->cval_alloc)
 147     {
 148         sfree(sb->cval);
 149         sb->cval_alloc = 0;
 150         sb->cval       = nullptr;
 151     }
 152     if (sb->sval_alloc)
 153     {
 154         int i;
 155
 156         for (i = 0; i < sb->sval_alloc; i++)
 157         {
 158             if (sb->sval[i])
 159             {
 160                 sfree(sb->sval[i]);
 161             }
 162         }
 163         sfree(sb->sval);
 164         sb->sval_alloc = 0;
 165         sb->sval       = nullptr;
 166     }
 167 }
 168
 169 /* allocate the appropriate amount of memory for the given type and nr */
 170 static void enxsubblock_alloc(t_enxsubblock *sb)
 171 {
 172     /* allocate the appropriate amount of memory */
 173     switch (sb->type)
 174     {
 175         case xdr_datatype_float:
 176             if (sb->nr > sb->fval_alloc)
 177             {
 178                 srenew(sb->fval, sb->nr);
 179                 sb->fval_alloc = sb->nr;
 180             }
 181             break;
 182         case xdr_datatype_double:
 183             if (sb->nr > sb->dval_alloc)
 184             {
 185                 srenew(sb->dval, sb->nr);
 186                 sb->dval_alloc = sb->nr;
 187             }
 188             break;
 189         case xdr_datatype_int:
 190             if (sb->nr > sb->ival_alloc)
 191             {
 192                 srenew(sb->ival, sb->nr);
 193                 sb->ival_alloc = sb->nr;
 194             }
 195             break;
 196         case xdr_datatype_int64:
 197             if (sb->nr > sb->lval_alloc)
 198             {
 199                 srenew(sb->lval, sb->nr);
 200                 sb->lval_alloc = sb->nr;
 201             }
 202             break;
 203         case xdr_datatype_char:
 204             if (sb->nr > sb->cval_alloc)
 205             {
 206                 srenew(sb->cval, sb->nr);
 207                 sb->cval_alloc = sb->nr;
 208             }
 209             break;
 210         case xdr_datatype_string:
 211             if (sb->nr > sb->sval_alloc)
 212             {
 213                 int i;
 214
 215                 srenew(sb->sval, sb->nr);
 216                 for (i = sb->sval_alloc; i < sb->nr; i++)
 217                 {
 218                     sb->sval[i] = nullptr;
 219                 }
 220                 sb->sval_alloc = sb->nr;
 221             }
 222             break;
 223         default:
 224             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(msg);
 336     }
 337     else
 338     {
 339         gmx_fatal(FARGS, "%s\n%s",
 340                   msg,
 341                   "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");
 342     }
 343 }
 344
 345 static void edr_strings(XDR *xdr, gmx_bool bRead, int file_version,
 346                         int n, gmx_enxnm_t **nms)
 347 {
 348     int          i;
 349     gmx_enxnm_t *nm;
 350
 351     if (*nms == nullptr)
 352     {
 353         snew(*nms, n);
 354     }
 355     for (i = 0; i < n; i++)
 356     {
 357         nm = &(*nms)[i];
 358         if (bRead)
 359         {
 360             if (nm->name)
 361             {
 362                 sfree(nm->name);
 363                 nm->name = nullptr;
 364             }
 365             if (nm->unit)
 366             {
 367                 sfree(nm->unit);
 368                 nm->unit = nullptr;
 369             }
 370         }
 371         if (!xdr_string(xdr, &(nm->name), STRLEN))
 372         {
 373             gmx_file("Cannot write energy names to file; maybe you are out of disk space?");
 374         }
 375         if (file_version >= 2)
 376         {
 377             if (!xdr_string(xdr, &(nm->unit), STRLEN))
 378             {
 379                 gmx_file("Cannot write energy names to file; maybe you are out of disk space?");
 380             }
 381         }
 382         else
 383         {
 384             nm->unit = gmx_strdup("kJ/mol");
 385         }
 386     }
 387 }
 388
 389 void do_enxnms(ener_file_t ef, int *nre, gmx_enxnm_t **nms)
 390 {
 391     int      magic = -55555;
 392     XDR     *xdr;
 393     gmx_bool bRead = gmx_fio_getread(ef->fio);
 394     int      file_version;
 395
 396     xdr = gmx_fio_getxdr(ef->fio);
 397
 398     if (!xdr_int(xdr, &magic))
 399     {
 400         if (!bRead)
 401         {
 402             gmx_file("Cannot write energy names to file; maybe you are out of disk space?");
 403         }
 404         *nre = 0;
 405         return;
 406     }
 407     if (magic > 0)
 408     {
 409         /* Assume this is an old edr format */
 410         file_version          = 1;
 411         *nre                  = magic;
 412         ef->eo.bOldFileOpen   = TRUE;
 413         ef->eo.bReadFirstStep = FALSE;
 414         srenew(ef->eo.ener_prev, *nre);
 415     }
 416     else
 417     {
 418         ef->eo.bOldFileOpen = FALSE;
 419
 420         if (magic != -55555)
 421         {
 422             gmx_fatal(FARGS, "Energy names magic number mismatch, this is not a GROMACS edr file");
 423         }
 424         file_version = enx_version;
 425         xdr_int(xdr, &file_version);
 426         if (file_version > enx_version)
 427         {
 428             gmx_fatal(FARGS, "reading tpx file (%s) version %d with version %d program", 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",
 435                 file_version, enx_version);
 436     }
 437
 438     edr_strings(xdr, bRead, file_version, *nre, nms);
 439 }
 440
 441 static gmx_bool do_eheader(ener_file_t ef, int *file_version, t_enxframe *fr,
 442                            int nre_test, gmx_bool *bWrongPrecision, gmx_bool *bOK)
 443 {
 444     int          magic = -7777777;
 445     real         first_real_to_check;
 446     int          b, zero = 0, dum = 0;
 447     gmx_bool     bRead      = gmx_fio_getread(ef->fio);
 448     int          ndisre     = 0;
 449     int          startb     = 0;
 450 #if !GMX_DOUBLE
 451     xdr_datatype dtreal = xdr_datatype_float;
 452 #else
 453     xdr_datatype dtreal = xdr_datatype_double;
 454 #endif
 455
 456     if (bWrongPrecision)
 457     {
 458         *bWrongPrecision = FALSE;
 459     }
 460
 461     *bOK = TRUE;
 462     /* The original energy frame started with a real,
 463      * so we have to use a real for compatibility.
 464      * This is VERY DIRTY code, since do_eheader can be called
 465      * with the wrong precision set and then we could read r > -1e10,
 466      * while actually the intention was r < -1e10.
 467      * When nre_test >= 0, do_eheader should therefore terminate
 468      * before the number of i/o calls starts depending on what has been read
 469      * (which is the case for for instance the block sizes for variable
 470      * number of blocks, where this number is read before).
 471      */
 472     first_real_to_check = -2e10;
 473     if (!gmx_fio_do_real(ef->fio, first_real_to_check))
 474     {
 475         return FALSE;
 476     }
 477     if (first_real_to_check > -1e10)
 478     {
 479         /* Assume we are reading an old format */
 480         *file_version = 1;
 481         fr->t         = first_real_to_check;
 482         if (!gmx_fio_do_int(ef->fio, dum))
 483         {
 484             *bOK = FALSE;
 485         }
 486         fr->step = dum;
 487     }
 488     else
 489     {
 490         if (!gmx_fio_do_int(ef->fio, magic))
 491         {
 492             *bOK = FALSE;
 493         }
 494         if (magic != -7777777)
 495         {
 496             enx_warning("Energy header magic number mismatch, this is not a GROMACS edr file");
 497             *bOK = FALSE;
 498             return FALSE;
 499         }
 500         *file_version = enx_version;
 501         // cppcheck-suppress redundantPointerOp
 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 = (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);
 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*(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*(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_t 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/(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=%" GMX_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     return 0;
1151 }
1152
1153
1154 void get_enx_state(const char *fn, real t, const gmx_groups_t *groups, t_inputrec *ir,
1155                    t_state *state)
1156 {
1157     /* Should match the names in mdebin.c */
1158     static const char *boxvel_nm[] = {
1159         "Box-Vel-XX", "Box-Vel-YY", "Box-Vel-ZZ",
1160         "Box-Vel-YX", "Box-Vel-ZX", "Box-Vel-ZY"
1161     };
1162
1163     static const char *baro_nm[] = {
1164         "Barostat"
1165     };
1166
1167
1168     int          ind0[] = { XX, YY, ZZ, YY, ZZ, ZZ };
1169     int          ind1[] = { XX, YY, ZZ, XX, XX, YY };
1170     int          nre, nfr, i, j, ni, npcoupl;
1171     char         buf[STRLEN];
1172     const char  *bufi;
1173     gmx_enxnm_t *enm = nullptr;
1174     t_enxframe  *fr;
1175     ener_file_t  in;
1176
1177     in = open_enx(fn, "r");
1178     do_enxnms(in, &nre, &enm);
1179     snew(fr, 1);
1180     nfr = 0;
1181     while ((nfr == 0 || fr->t != t) && do_enx(in, fr))
1182     {
1183         nfr++;
1184     }
1185     close_enx(in);
1186     fprintf(stderr, "\n");
1187
1188     if (nfr == 0 || fr->t != t)
1189     {
1190         gmx_fatal(FARGS, "Could not find frame with time %f in '%s'", t, fn);
1191     }
1192
1193     npcoupl = TRICLINIC(ir->compress) ? 6 : 3;
1194     if (ir->epc == epcPARRINELLORAHMAN)
1195     {
1196         clear_mat(state->boxv);
1197         for (i = 0; i < npcoupl; i++)
1198         {
1199             state->boxv[ind0[i]][ind1[i]] =
1200                 find_energy(boxvel_nm[i], nre, enm, fr);
1201         }
1202         fprintf(stderr, "\nREAD %d BOX VELOCITIES FROM %s\n\n", npcoupl, fn);
1203     }
1204
1205     if (ir->etc == etcNOSEHOOVER)
1206     {
1207         char cns[20];
1208
1209         cns[0] = '\0';
1210
1211         for (i = 0; i < state->ngtc; i++)
1212         {
1213             ni   = groups->grps[egcTC].nm_ind[i];
1214             bufi = *(groups->grpname[ni]);
1215             for (j = 0; (j < state->nhchainlength); j++)
1216             {
1217                 if (inputrecNvtTrotter(ir))
1218                 {
1219                     sprintf(cns, "-%d", j);
1220                 }
1221                 sprintf(buf, "Xi%s-%s", cns, bufi);
1222                 state->nosehoover_xi[i] = find_energy(buf, nre, enm, fr);
1223                 sprintf(buf, "vXi%s-%s", cns, bufi);
1224                 state->nosehoover_vxi[i] = find_energy(buf, nre, enm, fr);
1225             }
1226
1227         }
1228         fprintf(stderr, "\nREAD %d NOSE-HOOVER Xi chains FROM %s\n\n", state->ngtc, fn);
1229
1230         if (inputrecNptTrotter(ir) || inputrecNphTrotter(ir))
1231         {
1232             for (i = 0; i < state->nnhpres; i++)
1233             {
1234                 bufi = baro_nm[0]; /* All barostat DOF's together for now */
1235                 for (j = 0; (j < state->nhchainlength); j++)
1236                 {
1237                     sprintf(buf, "Xi-%d-%s", j, bufi);
1238                     state->nhpres_xi[i] = find_energy(buf, nre, enm, fr);
1239                     sprintf(buf, "vXi-%d-%s", j, bufi);
1240                     state->nhpres_vxi[i] = find_energy(buf, nre, enm, fr);
1241                 }
1242             }
1243             fprintf(stderr, "\nREAD %d NOSE-HOOVER BAROSTAT Xi chains FROM %s\n\n", state->nnhpres, fn);
1244         }
1245     }
1246
1247     free_enxnms(nre, enm);
1248     free_enxframe(fr);
1249     sfree(fr);
1250 }
1251
1252 static real ener_tensor_diag(int n, int *ind1, int *ind2,
1253                              gmx_enxnm_t *enm1,
1254                              int *tensi, int i,
1255                              t_energy e1[], t_energy e2[])
1256 {
1257     int    d1, d2;
1258     int    j;
1259     real   prod1, prod2;
1260     int    nfound;
1261     size_t len;
1262
1263     d1 = tensi[i]/DIM;
1264     d2 = tensi[i] - d1*DIM;
1265
1266     /* Find the diagonal elements d1 and d2 */
1267     len    = std::strlen(enm1[ind1[i]].name);
1268     prod1  = 1;
1269     prod2  = 1;
1270     nfound = 0;
1271     for (j = 0; j < n; j++)
1272     {
1273         if (tensi[j] >= 0 &&
1274             std::strlen(enm1[ind1[j]].name) == len &&
1275             std::strncmp(enm1[ind1[i]].name, enm1[ind1[j]].name, len-2) == 0 &&
1276             (tensi[j] == d1*DIM+d1 || tensi[j] == d2*DIM+d2))
1277         {
1278             prod1 *= fabs(e1[ind1[j]].e);
1279             prod2 *= fabs(e2[ind2[j]].e);
1280             nfound++;
1281         }
1282     }
1283
1284     if (nfound == 2)
1285     {
1286         return 0.5*(std::sqrt(prod1) + std::sqrt(prod2));
1287     }
1288     else
1289     {
1290         return 0;
1291     }
1292 }
1293
1294 static gmx_bool enernm_equal(const char *nm1, const char *nm2)
1295 {
1296     int len1, len2;
1297
1298     len1 = std::strlen(nm1);
1299     len2 = std::strlen(nm2);
1300
1301     /* Remove " (bar)" at the end of a name */
1302     if (len1 > 6 && std::strcmp(nm1+len1-6, " (bar)") == 0)
1303     {
1304         len1 -= 6;
1305     }
1306     if (len2 > 6 && std::strcmp(nm2+len2-6, " (bar)") == 0)
1307     {
1308         len2 -= 6;
1309     }
1310
1311     return (len1 == len2 && gmx_strncasecmp(nm1, nm2, len1) == 0);
1312 }
1313
1314 static void cmp_energies(FILE *fp, int step1, int step2,
1315                          t_energy e1[], t_energy e2[],
1316                          gmx_enxnm_t *enm1,
1317                          real ftol, real abstol,
1318                          int nre, int *ind1, int *ind2, int maxener)
1319 {
1320     int   i, ii;
1321     int  *tensi, len, d1, d2;
1322     real  ftol_i, abstol_i;
1323
1324     snew(tensi, maxener);
1325     /* Check for tensor elements ending on "-XX", "-XY", ... , "-ZZ" */
1326     for (i = 0; (i < maxener); i++)
1327     {
1328         ii       = ind1[i];
1329         tensi[i] = -1;
1330         len      = std::strlen(enm1[ii].name);
1331         if (len > 3 && enm1[ii].name[len-3] == '-')
1332         {
1333             d1 = enm1[ii].name[len-2] - 'X';
1334             d2 = enm1[ii].name[len-1] - 'X';
1335             if (d1 >= 0 && d1 < DIM &&
1336                 d2 >= 0 && d2 < DIM)
1337             {
1338                 tensi[i] = d1*DIM + d2;
1339             }
1340         }
1341     }
1342
1343     for (i = 0; (i < maxener); i++)
1344     {
1345         /* Check if this is an off-diagonal tensor element */
1346         if (tensi[i] >= 0 && tensi[i] != 0 && tensi[i] != 4 && tensi[i] != 8)
1347         {
1348             /* Turn on the relative tolerance check (4 is maximum relative diff.) */
1349             ftol_i = 5;
1350             /* Do the relative tolerance through an absolute tolerance times
1351              * the size of diagonal components of the tensor.
1352              */
1353             abstol_i = ftol*ener_tensor_diag(nre, ind1, ind2, enm1, tensi, i, e1, e2);
1354             if (debug)
1355             {
1356                 fprintf(debug, "tensor '%s' val %f diag %f\n",
1357                         enm1[i].name, e1[i].e, abstol_i/ftol);
1358             }
1359             if (abstol_i > 0)
1360             {
1361                 /* We found a diagonal, we need to check with the minimum tolerance */
1362                 abstol_i = std::min(abstol_i, abstol);
1363             }
1364             else
1365             {
1366                 /* We did not find a diagonal, ignore the relative tolerance check */
1367                 abstol_i = abstol;
1368             }
1369         }
1370         else
1371         {
1372             ftol_i   = ftol;
1373             abstol_i = abstol;
1374         }
1375         if (!equal_real(e1[ind1[i]].e, e2[ind2[i]].e, ftol_i, abstol_i))
1376         {
1377             fprintf(fp, "%-15s  step %3d:  %12g,  step %3d: %12g\n",
1378                     enm1[ind1[i]].name,
1379                     step1, e1[ind1[i]].e,
1380                     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, (int)s1->type, (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,
1446                                               s1->fval[k], s2->fval[k],
1447                                               ftol, abstol);
1448                                 }
1449                                 break;
1450                             case xdr_datatype_double:
1451                                 for (k = 0; k < s1->nr; k++)
1452                                 {
1453                                     cmp_double(stdout, buf, i,
1454                                                s1->dval[k], s2->dval[k],
1455                                                ftol, abstol);
1456                                 }
1457                                 break;
1458                             case xdr_datatype_int:
1459                                 for (k = 0; k < s1->nr; k++)
1460                                 {
1461                                     cmp_int(stdout, buf, i,
1462                                             s1->ival[k], s2->ival[k]);
1463                                 }
1464                                 break;
1465                             case xdr_datatype_int64:
1466                                 for (k = 0; k < s1->nr; k++)
1467                                 {
1468                                     cmp_int64(stdout, buf,
1469                                               s1->lval[k], s2->lval[k]);
1470                                 }
1471                                 break;
1472                             case xdr_datatype_char:
1473                                 for (k = 0; k < s1->nr; k++)
1474                                 {
1475                                     cmp_uc(stdout, buf, i,
1476                                            s1->cval[k], s2->cval[k]);
1477                                 }
1478                                 break;
1479                             case xdr_datatype_string:
1480                                 for (k = 0; k < s1->nr; k++)
1481                                 {
1482                                     cmp_str(stdout, buf, i,
1483                                             s1->sval[k], s2->sval[k]);
1484                                 }
1485                                 break;
1486                             default:
1487                                 gmx_incons("Unknown data type!!");
1488                         }
1489                     }
1490                 }
1491             }
1492         }
1493     }
1494 }
1495
1496 void comp_enx(const char *fn1, const char *fn2, real ftol, real abstol, const char *lastener)
1497 {
1498     int            nre, nre1, nre2;
1499     ener_file_t    in1, in2;
1500     int            i, j, maxener, *ind1, *ind2, *have;
1501     gmx_enxnm_t   *enm1 = nullptr, *enm2 = nullptr;
1502     t_enxframe    *fr1, *fr2;
1503     gmx_bool       b1, b2;
1504
1505     fprintf(stdout, "comparing energy file %s and %s\n\n", fn1, fn2);
1506
1507     in1 = open_enx(fn1, "r");
1508     in2 = open_enx(fn2, "r");
1509     do_enxnms(in1, &nre1, &enm1);
1510     do_enxnms(in2, &nre2, &enm2);
1511     if (nre1 != nre2)
1512     {
1513         fprintf(stdout, "There are %d and %d terms in the energy files\n\n",
1514                 nre1, nre2);
1515     }
1516     else
1517     {
1518         fprintf(stdout, "There are %d terms in the energy files\n\n", nre1);
1519     }
1520
1521     snew(ind1, nre1);
1522     snew(ind2, nre2);
1523     snew(have, nre2);
1524     nre = 0;
1525     for (i = 0; i < nre1; i++)
1526     {
1527         for (j = 0; j < nre2; j++)
1528         {
1529             if (enernm_equal(enm1[i].name, enm2[j].name))
1530             {
1531                 ind1[nre] = i;
1532                 ind2[nre] = j;
1533                 have[j]   = 1;
1534                 nre++;
1535                 break;
1536             }
1537         }
1538         if (nre == 0 || ind1[nre-1] != i)
1539         {
1540             cmp_str(stdout, "enm", i, enm1[i].name, "-");
1541         }
1542     }
1543     for (i = 0; i < nre2; i++)
1544     {
1545         if (have[i] == 0)
1546         {
1547             cmp_str(stdout, "enm", i, "-", enm2[i].name);
1548         }
1549     }
1550
1551     maxener = nre;
1552     for (i = 0; i < nre; i++)
1553     {
1554         if ((lastener != nullptr) && (std::strstr(enm1[i].name, lastener) != nullptr))
1555         {
1556             maxener = i+1;
1557             break;
1558         }
1559     }
1560
1561     fprintf(stdout, "There are %d terms to compare in the energy files\n\n",
1562             maxener);
1563
1564     for (i = 0; i < maxener; i++)
1565     {
1566         cmp_str(stdout, "unit", i, enm1[ind1[i]].unit, enm2[ind2[i]].unit);
1567     }
1568
1569     snew(fr1, 1);
1570     snew(fr2, 1);
1571     do
1572     {
1573         b1 = do_enx(in1, fr1);
1574         b2 = do_enx(in2, fr2);
1575         if (b1 && !b2)
1576         {
1577             fprintf(stdout, "\nEnd of file on %s but not on %s\n", fn2, fn1);
1578         }
1579         else if (!b1 && b2)
1580         {
1581             fprintf(stdout, "\nEnd of file on %s but not on %s\n", fn1, fn2);
1582         }
1583         else if (!b1 && !b2)
1584         {
1585             fprintf(stdout, "\nFiles read successfully\n");
1586         }
1587         else
1588         {
1589             cmp_real(stdout, "t", -1, fr1->t, fr2->t, ftol, abstol);
1590             cmp_int(stdout, "step", -1, fr1->step, fr2->step);
1591             /* We don't want to print the nre mismatch for every frame */
1592             /* cmp_int(stdout,"nre",-1,fr1->nre,fr2->nre); */
1593             if ((fr1->nre >= nre) && (fr2->nre >= nre))
1594             {
1595                 cmp_energies(stdout, fr1->step, fr1->step, fr1->ener, fr2->ener,
1596                              enm1, ftol, abstol, nre, ind1, ind2, maxener);
1597             }
1598             /*cmp_disres(fr1,fr2,ftol,abstol);*/
1599             cmp_eblocks(fr1, fr2, ftol, abstol);
1600         }
1601     }
1602     while (b1 && b2);
1603
1604     close_enx(in1);
1605     close_enx(in2);
1606
1607     free_enxframe(fr2);
1608     sfree(fr2);
1609     free_enxframe(fr1);
1610     sfree(fr1);
1611 }