src/gromacs/fileio/enxio.cpp

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