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