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(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 = (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*(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*(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_t 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/(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=%" GMX_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     return 0;
1152 }
1153
1154
1155 void get_enx_state(const char *fn, real t, const gmx_groups_t *groups, t_inputrec *ir,
1156                    t_state *state)
1157 {
1158     /* Should match the names in mdebin.c */
1159     static const char *boxvel_nm[] = {
1160         "Box-Vel-XX", "Box-Vel-YY", "Box-Vel-ZZ",
1161         "Box-Vel-YX", "Box-Vel-ZX", "Box-Vel-ZY"
1162     };
1163
1164     static const char *baro_nm[] = {
1165         "Barostat"
1166     };
1167
1168
1169     int          ind0[] = { XX, YY, ZZ, YY, ZZ, ZZ };
1170     int          ind1[] = { XX, YY, ZZ, XX, XX, YY };
1171     int          nre, nfr, i, j, ni, npcoupl;
1172     char         buf[STRLEN];
1173     const char  *bufi;
1174     gmx_enxnm_t *enm = nullptr;
1175     t_enxframe  *fr;
1176     ener_file_t  in;
1177
1178     in = open_enx(fn, "r");
1179     do_enxnms(in, &nre, &enm);
1180     snew(fr, 1);
1181     nfr = 0;
1182     while ((nfr == 0 || fr->t != t) && do_enx(in, fr))
1183     {
1184         nfr++;
1185     }
1186     close_enx(in);
1187     fprintf(stderr, "\n");
1188
1189     if (nfr == 0 || fr->t != t)
1190     {
1191         gmx_fatal(FARGS, "Could not find frame with time %f in '%s'", t, fn);
1192     }
1193
1194     npcoupl = TRICLINIC(ir->compress) ? 6 : 3;
1195     if (ir->epc == epcPARRINELLORAHMAN)
1196     {
1197         clear_mat(state->boxv);
1198         for (i = 0; i < npcoupl; i++)
1199         {
1200             state->boxv[ind0[i]][ind1[i]] =
1201                 find_energy(boxvel_nm[i], nre, enm, fr);
1202         }
1203         fprintf(stderr, "\nREAD %d BOX VELOCITIES FROM %s\n\n", npcoupl, fn);
1204     }
1205
1206     if (ir->etc == etcNOSEHOOVER)
1207     {
1208         char cns[20];
1209
1210         cns[0] = '\0';
1211
1212         for (i = 0; i < state->ngtc; i++)
1213         {
1214             ni   = groups->grps[egcTC].nm_ind[i];
1215             bufi = *(groups->grpname[ni]);
1216             for (j = 0; (j < state->nhchainlength); j++)
1217             {
1218                 if (inputrecNvtTrotter(ir))
1219                 {
1220                     sprintf(cns, "-%d", j);
1221                 }
1222                 sprintf(buf, "Xi%s-%s", cns, bufi);
1223                 state->nosehoover_xi[i] = find_energy(buf, nre, enm, fr);
1224                 sprintf(buf, "vXi%s-%s", cns, bufi);
1225                 state->nosehoover_vxi[i] = find_energy(buf, nre, enm, fr);
1226             }
1227
1228         }
1229         fprintf(stderr, "\nREAD %d NOSE-HOOVER Xi chains FROM %s\n\n", state->ngtc, fn);
1230
1231         if (inputrecNptTrotter(ir) || inputrecNphTrotter(ir))
1232         {
1233             for (i = 0; i < state->nnhpres; i++)
1234             {
1235                 bufi = baro_nm[0]; /* All barostat DOF's together for now */
1236                 for (j = 0; (j < state->nhchainlength); j++)
1237                 {
1238                     sprintf(buf, "Xi-%d-%s", j, bufi);
1239                     state->nhpres_xi[i] = find_energy(buf, nre, enm, fr);
1240                     sprintf(buf, "vXi-%d-%s", j, bufi);
1241                     state->nhpres_vxi[i] = find_energy(buf, nre, enm, fr);
1242                 }
1243             }
1244             fprintf(stderr, "\nREAD %d NOSE-HOOVER BAROSTAT Xi chains FROM %s\n\n", state->nnhpres, fn);
1245         }
1246     }
1247
1248     free_enxnms(nre, enm);
1249     free_enxframe(fr);
1250     sfree(fr);
1251 }
1252
1253 static real ener_tensor_diag(int n, int *ind1, int *ind2,
1254                              gmx_enxnm_t *enm1,
1255                              int *tensi, int i,
1256                              t_energy e1[], t_energy e2[])
1257 {
1258     int    d1, d2;
1259     int    j;
1260     real   prod1, prod2;
1261     int    nfound;
1262     size_t len;
1263
1264     d1 = tensi[i]/DIM;
1265     d2 = tensi[i] - d1*DIM;
1266
1267     /* Find the diagonal elements d1 and d2 */
1268     len    = std::strlen(enm1[ind1[i]].name);
1269     prod1  = 1;
1270     prod2  = 1;
1271     nfound = 0;
1272     for (j = 0; j < n; j++)
1273     {
1274         if (tensi[j] >= 0 &&
1275             std::strlen(enm1[ind1[j]].name) == len &&
1276             std::strncmp(enm1[ind1[i]].name, enm1[ind1[j]].name, len-2) == 0 &&
1277             (tensi[j] == d1*DIM+d1 || tensi[j] == d2*DIM+d2))
1278         {
1279             prod1 *= fabs(e1[ind1[j]].e);
1280             prod2 *= fabs(e2[ind2[j]].e);
1281             nfound++;
1282         }
1283     }
1284
1285     if (nfound == 2)
1286     {
1287         return 0.5*(std::sqrt(prod1) + std::sqrt(prod2));
1288     }
1289     else
1290     {
1291         return 0;
1292     }
1293 }
1294
1295 static gmx_bool enernm_equal(const char *nm1, const char *nm2)
1296 {
1297     int len1, len2;
1298
1299     len1 = std::strlen(nm1);
1300     len2 = std::strlen(nm2);
1301
1302     /* Remove " (bar)" at the end of a name */
1303     if (len1 > 6 && std::strcmp(nm1+len1-6, " (bar)") == 0)
1304     {
1305         len1 -= 6;
1306     }
1307     if (len2 > 6 && std::strcmp(nm2+len2-6, " (bar)") == 0)
1308     {
1309         len2 -= 6;
1310     }
1311
1312     return (len1 == len2 && gmx_strncasecmp(nm1, nm2, len1) == 0);
1313 }
1314
1315 static void cmp_energies(FILE *fp, int step1, int step2,
1316                          t_energy e1[], t_energy e2[],
1317                          gmx_enxnm_t *enm1,
1318                          real ftol, real abstol,
1319                          int nre, int *ind1, int *ind2, int maxener)
1320 {
1321     int   i, ii;
1322     int  *tensi, len, d1, d2;
1323     real  ftol_i, abstol_i;
1324
1325     snew(tensi, maxener);
1326     /* Check for tensor elements ending on "-XX", "-XY", ... , "-ZZ" */
1327     for (i = 0; (i < maxener); i++)
1328     {
1329         ii       = ind1[i];
1330         tensi[i] = -1;
1331         len      = std::strlen(enm1[ii].name);
1332         if (len > 3 && enm1[ii].name[len-3] == '-')
1333         {
1334             d1 = enm1[ii].name[len-2] - 'X';
1335             d2 = enm1[ii].name[len-1] - 'X';
1336             if (d1 >= 0 && d1 < DIM &&
1337                 d2 >= 0 && d2 < DIM)
1338             {
1339                 tensi[i] = d1*DIM + d2;
1340             }
1341         }
1342     }
1343
1344     for (i = 0; (i < maxener); i++)
1345     {
1346         /* Check if this is an off-diagonal tensor element */
1347         if (tensi[i] >= 0 && tensi[i] != 0 && tensi[i] != 4 && tensi[i] != 8)
1348         {
1349             /* Turn on the relative tolerance check (4 is maximum relative diff.) */
1350             ftol_i = 5;
1351             /* Do the relative tolerance through an absolute tolerance times
1352              * the size of diagonal components of the tensor.
1353              */
1354             abstol_i = ftol*ener_tensor_diag(nre, ind1, ind2, enm1, tensi, i, e1, e2);
1355             if (debug)
1356             {
1357                 fprintf(debug, "tensor '%s' val %f diag %f\n",
1358                         enm1[i].name, e1[i].e, abstol_i/ftol);
1359             }
1360             if (abstol_i > 0)
1361             {
1362                 /* We found a diagonal, we need to check with the minimum tolerance */
1363                 abstol_i = std::min(abstol_i, abstol);
1364             }
1365             else
1366             {
1367                 /* We did not find a diagonal, ignore the relative tolerance check */
1368                 abstol_i = abstol;
1369             }
1370         }
1371         else
1372         {
1373             ftol_i   = ftol;
1374             abstol_i = abstol;
1375         }
1376         if (!equal_real(e1[ind1[i]].e, e2[ind2[i]].e, ftol_i, abstol_i))
1377         {
1378             fprintf(fp, "%-15s  step %3d:  %12g,  step %3d: %12g\n",
1379                     enm1[ind1[i]].name,
1380                     step1, e1[ind1[i]].e,
1381                     step2, e2[ind2[i]].e);
1382         }
1383     }
1384
1385     sfree(tensi);
1386 }
1387
1388 #if 0
1389 static void cmp_disres(t_enxframe *fr1, t_enxframe *fr2, real ftol, real abstol)
1390 {
1391     int  i;
1392     char bav[64], bt[64], bs[22];
1393
1394     cmp_int(stdout, "ndisre", -1, fr1->ndisre, fr2->ndisre);
1395     if ((fr1->ndisre == fr2->ndisre) && (fr1->ndisre > 0))
1396     {
1397         sprintf(bav, "step %s: disre rav", gmx_step_str(fr1->step, bs));
1398         sprintf(bt, "step %s: disre  rt", gmx_step_str(fr1->step, bs));
1399         for (i = 0; (i < fr1->ndisre); i++)
1400         {
1401             cmp_real(stdout, bav, i, fr1->disre_rm3tav[i], fr2->disre_rm3tav[i], ftol, abstol);
1402             cmp_real(stdout, bt, i, fr1->disre_rt[i], fr2->disre_rt[i], ftol, abstol);
1403         }
1404     }
1405 }
1406 #endif
1407
1408 static void cmp_eblocks(t_enxframe *fr1, t_enxframe *fr2, real ftol, real abstol)
1409 {
1410     int  i, j, k;
1411     char buf[64], bs[22];
1412
1413     cmp_int(stdout, "nblock", -1, fr1->nblock, fr2->nblock);
1414     if ((fr1->nblock == fr2->nblock) && (fr1->nblock > 0))
1415     {
1416         for (j = 0; (j < fr1->nblock); j++)
1417         {
1418             t_enxblock *b1, *b2; /* convenience vars */
1419
1420             b1 = &(fr1->block[j]);
1421             b2 = &(fr2->block[j]);
1422
1423             sprintf(buf, "step %s: block[%d]", gmx_step_str(fr1->step, bs), j);
1424             cmp_int(stdout, buf, -1, b1->nsub, b2->nsub);
1425             cmp_int(stdout, buf, -1, b1->id, b2->id);
1426
1427             if ( (b1->nsub == b2->nsub) && (b1->id == b2->id) )
1428             {
1429                 for (i = 0; i < b1->nsub; i++)
1430                 {
1431                     t_enxsubblock *s1, *s2;
1432
1433                     s1 = &(b1->sub[i]);
1434                     s2 = &(b2->sub[i]);
1435
1436                     cmp_int(stdout, buf, -1, (int)s1->type, (int)s2->type);
1437                     cmp_int64(stdout, buf, s1->nr, s2->nr);
1438
1439                     if ((s1->type == s2->type) && (s1->nr == s2->nr))
1440                     {
1441                         switch (s1->type)
1442                         {
1443                             case xdr_datatype_float:
1444                                 for (k = 0; k < s1->nr; k++)
1445                                 {
1446                                     cmp_float(stdout, buf, i,
1447                                               s1->fval[k], s2->fval[k],
1448                                               ftol, abstol);
1449                                 }
1450                                 break;
1451                             case xdr_datatype_double:
1452                                 for (k = 0; k < s1->nr; k++)
1453                                 {
1454                                     cmp_double(stdout, buf, i,
1455                                                s1->dval[k], s2->dval[k],
1456                                                ftol, abstol);
1457                                 }
1458                                 break;
1459                             case xdr_datatype_int:
1460                                 for (k = 0; k < s1->nr; k++)
1461                                 {
1462                                     cmp_int(stdout, buf, i,
1463                                             s1->ival[k], s2->ival[k]);
1464                                 }
1465                                 break;
1466                             case xdr_datatype_int64:
1467                                 for (k = 0; k < s1->nr; k++)
1468                                 {
1469                                     cmp_int64(stdout, buf,
1470                                               s1->lval[k], s2->lval[k]);
1471                                 }
1472                                 break;
1473                             case xdr_datatype_char:
1474                                 for (k = 0; k < s1->nr; k++)
1475                                 {
1476                                     cmp_uc(stdout, buf, i,
1477                                            s1->cval[k], s2->cval[k]);
1478                                 }
1479                                 break;
1480                             case xdr_datatype_string:
1481                                 for (k = 0; k < s1->nr; k++)
1482                                 {
1483                                     cmp_str(stdout, buf, i,
1484                                             s1->sval[k], s2->sval[k]);
1485                                 }
1486                                 break;
1487                             default:
1488                                 gmx_incons("Unknown data type!!");
1489                         }
1490                     }
1491                 }
1492             }
1493         }
1494     }
1495 }
1496
1497 void comp_enx(const char *fn1, const char *fn2, real ftol, real abstol, const char *lastener)
1498 {
1499     int            nre, nre1, nre2;
1500     ener_file_t    in1, in2;
1501     int            i, j, maxener, *ind1, *ind2, *have;
1502     gmx_enxnm_t   *enm1 = nullptr, *enm2 = nullptr;
1503     t_enxframe    *fr1, *fr2;
1504     gmx_bool       b1, b2;
1505
1506     fprintf(stdout, "comparing energy file %s and %s\n\n", fn1, fn2);
1507
1508     in1 = open_enx(fn1, "r");
1509     in2 = open_enx(fn2, "r");
1510     do_enxnms(in1, &nre1, &enm1);
1511     do_enxnms(in2, &nre2, &enm2);
1512     if (nre1 != nre2)
1513     {
1514         fprintf(stdout, "There are %d and %d terms in the energy files\n\n",
1515                 nre1, nre2);
1516     }
1517     else
1518     {
1519         fprintf(stdout, "There are %d terms in the energy files\n\n", nre1);
1520     }
1521
1522     snew(ind1, nre1);
1523     snew(ind2, nre2);
1524     snew(have, nre2);
1525     nre = 0;
1526     for (i = 0; i < nre1; i++)
1527     {
1528         for (j = 0; j < nre2; j++)
1529         {
1530             if (enernm_equal(enm1[i].name, enm2[j].name))
1531             {
1532                 ind1[nre] = i;
1533                 ind2[nre] = j;
1534                 have[j]   = 1;
1535                 nre++;
1536                 break;
1537             }
1538         }
1539         if (nre == 0 || ind1[nre-1] != i)
1540         {
1541             cmp_str(stdout, "enm", i, enm1[i].name, "-");
1542         }
1543     }
1544     for (i = 0; i < nre2; i++)
1545     {
1546         if (have[i] == 0)
1547         {
1548             cmp_str(stdout, "enm", i, "-", enm2[i].name);
1549         }
1550     }
1551
1552     maxener = nre;
1553     for (i = 0; i < nre; i++)
1554     {
1555         if ((lastener != nullptr) && (std::strstr(enm1[i].name, lastener) != nullptr))
1556         {
1557             maxener = i+1;
1558             break;
1559         }
1560     }
1561
1562     fprintf(stdout, "There are %d terms to compare in the energy files\n\n",
1563             maxener);
1564
1565     for (i = 0; i < maxener; i++)
1566     {
1567         cmp_str(stdout, "unit", i, enm1[ind1[i]].unit, enm2[ind2[i]].unit);
1568     }
1569
1570     snew(fr1, 1);
1571     snew(fr2, 1);
1572     do
1573     {
1574         b1 = do_enx(in1, fr1);
1575         b2 = do_enx(in2, fr2);
1576         if (b1 && !b2)
1577         {
1578             fprintf(stdout, "\nEnd of file on %s but not on %s\n", fn2, fn1);
1579         }
1580         else if (!b1 && b2)
1581         {
1582             fprintf(stdout, "\nEnd of file on %s but not on %s\n", fn1, fn2);
1583         }
1584         else if (!b1 && !b2)
1585         {
1586             fprintf(stdout, "\nFiles read successfully\n");
1587         }
1588         else
1589         {
1590             cmp_real(stdout, "t", -1, fr1->t, fr2->t, ftol, abstol);
1591             cmp_int(stdout, "step", -1, fr1->step, fr2->step);
1592             /* We don't want to print the nre mismatch for every frame */
1593             /* cmp_int(stdout,"nre",-1,fr1->nre,fr2->nre); */
1594             if ((fr1->nre >= nre) && (fr2->nre >= nre))
1595             {
1596                 cmp_energies(stdout, fr1->step, fr1->step, fr1->ener, fr2->ener,
1597                              enm1, ftol, abstol, nre, ind1, ind2, maxener);
1598             }
1599             /*cmp_disres(fr1,fr2,ftol,abstol);*/
1600             cmp_eblocks(fr1, fr2, ftol, abstol);
1601         }
1602     }
1603     while (b1 && b2);
1604
1605     close_enx(in1);
1606     close_enx(in2);
1607
1608     free_enxframe(fr2);
1609     sfree(fr2);
1610     free_enxframe(fr1);
1611     sfree(fr1);
1612 }