src/gromacs/gmxana/gmx_awh.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35
  36 /*! \internal \file
  37  *  \brief
  38  *  Tool for extracting AWH (Accelerated Weight Histogram) data from energy files.
  39  *
  40  *  \author Viveca Lindahl
  41  *  \author Berk Hess
  42  */
  43
  44 #include "gmxpre.h"
  45
  46 #include <cstdlib>
  47 #include <cstring>
  48
  49 #include <algorithm>
  50 #include <array>
  51 #include <memory>
  52 #include <string>
  53
  54 #include "gromacs/commandline/pargs.h"
  55 #include "gromacs/compat/make_unique.h"
  56 #include "gromacs/fileio/enxio.h"
  57 #include "gromacs/fileio/gmxfio.h"
  58 #include "gromacs/fileio/oenv.h"
  59 #include "gromacs/fileio/tpxio.h"
  60 #include "gromacs/fileio/trxio.h"
  61 #include "gromacs/fileio/xvgr.h"
  62 #include "gromacs/gmxana/gmx_ana.h"
  63 #include "gromacs/math/units.h"
  64 #include "gromacs/mdtypes/awh-params.h"
  65 #include "gromacs/mdtypes/inputrec.h"
  66 #include "gromacs/topology/mtop_util.h"
  67 #include "gromacs/topology/topology.h"
  68 #include "gromacs/trajectory/energyframe.h"
  69 #include "gromacs/utility/arraysize.h"
  70 #include "gromacs/utility/basedefinitions.h"
  71 #include "gromacs/utility/fatalerror.h"
  72 #include "gromacs/utility/futil.h"
  73 #include "gromacs/utility/gmxassert.h"
  74 #include "gromacs/utility/smalloc.h"
  75 #include "gromacs/utility/stringutil.h"
  76
  77 using gmx::AwhParams;
  78 using gmx::AwhBiasParams;
  79
  80 namespace
  81 {
  82
  83 //! Enum for choosing AWH graph output
  84 enum class AwhGraphSelection
  85 {
  86     Pmf, //!< Only the PMF
  87     All  //!< All possible AWH graphs
  88 };
  89
  90 //! Energy unit options
  91 enum class EnergyUnit
  92 {
  93     KJPerMol, //!< kJ/mol
  94     KT        //!< kT
  95 };
  96
  97 } // namespace
  98
  99 /*! \brief All meta-data that is shared for one output file type for one bias */
 100 class OutputFile
 101 {
 102     public:
 103         /*! \brief Constructor, Set the output base file name and title.
 104          *
 105          * Result is a valid object, but will produce empty output files.
 106          *
 107          * \param[in] filename   The name for output files, frame time, and possibly bias number, will be added per file/frame.
 108          * \param[in] baseTitle  The base title of the plot, the bias number might be added.
 109          * \param[in] numBias    The total number of AWH biases in the system.
 110          * \param[in] biasIndex  The index of this bias.
 111          */
 112         OutputFile(const std::string  filename,
 113                    const std::string  baseTitle,
 114                    int                numBias,
 115                    int                biasIndex);
 116
 117         /*! \brief Initializes the output file setup for the AWH output.
 118          *
 119          * \param[in] subBlockStart   Index of the first sub-block to write in the energy frame.
 120          * \param[in] numSubBlocks    The total number of sub-blocks in the framw.
 121          * \param[in] awhBiasParams   The AWH bias parameters.
 122          * \param[in] graphSelection  Selects which graphs to plot.
 123          * \param[in] energyUnit      Requested energy unit in output.
 124          * \param[in] kTValue         kB*T in kJ/mol.
 125          */
 126         void initializeAwhOutputFile(int                  subBlockStart,
 127                                      int                  numSubBlocks,
 128                                      const AwhBiasParams *awhBiasParams,
 129                                      AwhGraphSelection    graphSelection,
 130                                      EnergyUnit           energyUnit,
 131                                      real                 kTValue);
 132
 133         /*! \brief Initializes the output file setup for the fricion output.
 134          *
 135          * \param[in] subBlockStart   Index of the first sub-block to write in the energy frame.
 136          * \param[in] numSubBlocks    The total number of sub-blocks in the framw.
 137          * \param[in] awhBiasParams   The AWH bias parameters.
 138          * \param[in] energyUnit      Requested energy unit in output.
 139          * \param[in] kTValue         kB*T in kJ/mol.
 140          */
 141         void initializeFrictionOutputFile(int                  subBlockStart,
 142                                           int                  numSubBlocks,
 143                                           const AwhBiasParams *awhBiasParams,
 144                                           EnergyUnit           energyUnit,
 145                                           real                 kTValue);
 146
 147         /*! \brief Opens a single output file for a bias, prints title and legends.
 148          *
 149          * \param[in] time  The time for of frame to be written.
 150          * \param[in] oenv  The output environment.
 151          * \returns the file pointer.
 152          */
 153         FILE * openBiasOutputFile(double                  time,
 154                                   const gmx_output_env_t *oenv) const;
 155
 156         /*! \brief Writes data selected from \p block to file.
 157          *
 158          * \param[in] block          The energy block with the data to write.
 159          * \param[in] subBlockStart  The sub-block start index.
 160          * \param[in] fp             The file pointer.
 161          */
 162         void writeData(const t_enxblock &block,
 163                        int               subBlockStart,
 164                        FILE             *fp) const;
 165
 166     private:
 167         std::string              baseFilename_;       /**< Base of the output file name. */
 168         std::string              title_;              /**< Title for the graph. */
 169         int                      numDim_;             /**< Number of dimensions. */
 170         int                      firstGraphSubBlock_; /**< Index in the energy sub-blocks for the first graph. */
 171         int                      numGraph_;           /**< Number of actual graphs. */
 172         bool                     useKTForEnergy_;     /**< Whether to use kT instead of kJ/mol for energy. */
 173         std::vector<real>        scaleFactor_;        /**< Scale factors from energy file data to output for each of the numGraph quantities. */
 174
 175         std::vector<std::string> legend_;             /**< Legends for the output */
 176         std::string              xLabel_;             /**< Label for the x-axis. */
 177         std::string              yLabel_;             /**< Label for the y-axis. */
 178 };
 179
 180 /*! \brief All meta-data that is shared for all output files for one bias */
 181 struct BiasOutputSetup
 182 {
 183     //! Constructor.
 184     BiasOutputSetup(int                         subblockStart,
 185                     std::unique_ptr<OutputFile> awhOutputFile,
 186                     std::unique_ptr<OutputFile> frictionOutputFile) :
 187         subblockStart_(subblockStart),
 188         awhOutputFile_(std::move(awhOutputFile)),
 189         frictionOutputFile_(std::move(frictionOutputFile))
 190     {
 191     }
 192
 193     //! Return the AWH output file data.
 194     const OutputFile &awhOutputFile() const
 195     {
 196         GMX_RELEASE_ASSERT(awhOutputFile_ != nullptr, "awhOutputFile() called without initialized AWH output file");
 197
 198         return *awhOutputFile_.get();
 199     }
 200
 201     //! Return the a pointer to the friction output file data, can return nullptr
 202     const OutputFile *frictionOutputFile() const
 203     {
 204         return frictionOutputFile_.get();
 205     }
 206
 207     //! Return the starting subblock.
 208     int subblockStart() const
 209     {
 210         return subblockStart_;
 211     }
 212
 213     private:
 214         const int                   subblockStart_;      /**< The start index of the subblocks to read. */
 215         std::unique_ptr<OutputFile> awhOutputFile_;      /**< The standard AWH output file data. */
 216         std::unique_ptr<OutputFile> frictionOutputFile_; /**< The friction/metric tensor output file data */
 217 };
 218
 219 /*! \brief All options and meta-data needed for the AWH output */
 220 class AwhReader
 221 {
 222     public:
 223         //! Constructor
 224         AwhReader(const AwhParams  *awhParams,
 225                   int               numFileOptions,
 226                   const t_filenm   *filenames,
 227                   AwhGraphSelection awhGraphSelection,
 228                   EnergyUnit        energyUnit,
 229                   real              kT,
 230                   const t_enxblock *block);
 231
 232         //! Extract and print AWH data for an AWH block of one energy frame
 233         void processAwhFrame(const t_enxblock       &block,
 234                              double                  time,
 235                              const gmx_output_env_t *oenv) const;
 236
 237     private:
 238         std::vector<BiasOutputSetup> biasOutputSetups_; /**< The output setups, one for each AWH bias. */
 239     public:
 240         const real                   kT_;               /**< kB*T in kJ/mol. */
 241 };
 242
 243 namespace
 244 {
 245
 246 //! Enum for selecting output file type.
 247 enum class OutputFileType
 248 {
 249     Awh,      //!< AWH, all data except friction tensor.
 250     Friction  //!< The full friction tensor.
 251 };
 252
 253 /*! \brief The maximum number of observables per subblock, not including the full friction tensor.
 254  *
 255  * This number, as well as the legends in make_legend() should match
 256  * the output that mdrun writes. It would be better to define these
 257  * values in a single location.
 258  */
 259 static constexpr int maxAwhGraphs = 6;
 260
 261 /*! \brief Constructs a legend for a standard awh output file */
 262 std::vector<std::string>makeLegend(const AwhBiasParams *awhBiasParams,
 263                                    OutputFileType       outputFileType,
 264                                    size_t               numLegend)
 265 {
 266     const std::array<std::string, maxAwhGraphs> legendBase =
 267     {
 268         { "PMF",
 269           "Coord bias",
 270           "Coord distr",
 271           "Ref value distr",
 272           "Target ref value distr",
 273           "Friction metric" }
 274     };
 275
 276     std::vector<std::string>                    legend;
 277     /* Give legends to dimensions higher than the first */
 278     for (int d = 1; d < awhBiasParams->ndim; d++)
 279     {
 280         legend.push_back(gmx::formatString("dim%d", d));
 281     }
 282
 283     switch (outputFileType)
 284     {
 285         case OutputFileType::Awh:
 286         {
 287             /* Add as many legends as possible from the "base" legend list */
 288             size_t legendBaseIndex = 0;
 289             while (legend.size() < numLegend && legendBaseIndex < legendBase.size())
 290             {
 291                 legend.push_back(legendBase[legendBaseIndex]);
 292                 legendBaseIndex++;
 293             }
 294         }
 295         break;
 296         case OutputFileType::Friction:
 297             for (int i0 = 0; i0 < awhBiasParams->ndim; i0++)
 298             {
 299                 for (int i1 = 0; i1 <= i0; i1++)
 300                 {
 301                     legend.push_back(gmx::formatString("%d,%d", i0, i1));
 302                 }
 303             }
 304             break;
 305     }
 306
 307     GMX_RELEASE_ASSERT(legend.size() == numLegend, "The number of legends requested for printing and the number generated should be the same");
 308
 309     return legend;
 310 }
 311
 312 } // namespace
 313
 314 OutputFile::OutputFile(const std::string  filename,
 315                        const std::string  baseTitle,
 316                        int                numBias,
 317                        int                biasIndex) :
 318     numDim_(0),
 319     firstGraphSubBlock_(0),
 320     numGraph_(0),
 321     useKTForEnergy_(0),
 322     scaleFactor_(),
 323     legend_(),
 324     xLabel_(),
 325     yLabel_()
 326 {
 327     // cppcheck-suppress useInitializationList
 328     baseFilename_ = filename.substr(0, filename.find('.'));
 329     title_        = baseTitle;
 330     if (numBias > 1)
 331     {
 332         baseFilename_ += gmx::formatString("%d", biasIndex + 1);
 333         title_        += gmx::formatString(" %d", biasIndex + 1);
 334     }
 335 }
 336
 337 void OutputFile::initializeAwhOutputFile(int                  subblockStart,
 338                                          int                  numSubBlocks,
 339                                          const AwhBiasParams *awhBiasParams,
 340                                          AwhGraphSelection    graphSelection,
 341                                          EnergyUnit           energyUnit,
 342                                          real                 kTValue)
 343 {
 344     /* The first subblock with actual graph y-values is index 1 + ndim */
 345     numDim_             = awhBiasParams->ndim;
 346     firstGraphSubBlock_ = subblockStart + 1 + numDim_;
 347     if (graphSelection == AwhGraphSelection::All)
 348     {
 349         /* There are one metadata and ndim coordinate blocks */
 350         numGraph_       = std::min(numSubBlocks - 1 - numDim_,
 351                                    maxAwhGraphs);
 352     }
 353     else
 354     {
 355         numGraph_       = 1;
 356     }
 357     useKTForEnergy_     = (energyUnit == EnergyUnit::KT);
 358     scaleFactor_.resize(numGraph_, 1);
 359     if (!useKTForEnergy_)
 360     {
 361         /* The first two graphs are in units of energy, multiply by kT */
 362         std::fill(scaleFactor_.begin(), scaleFactor_.begin() + std::min(2, numGraph_), kTValue);
 363     }
 364     int numLegend = numDim_ - 1 + numGraph_;
 365     legend_       = makeLegend(awhBiasParams, OutputFileType::Awh, numLegend);
 366     /* We could have both length and angle coordinates in a single bias */
 367     xLabel_       = "(nm or deg)";
 368     yLabel_       = useKTForEnergy_ ? "(k\\sB\\NT)" : "(kJ/mol)";
 369     if (graphSelection == AwhGraphSelection::All)
 370     {
 371         yLabel_ += gmx::formatString(", (nm\\S-%d\\N or rad\\S-%d\\N), (-)",
 372                                      awhBiasParams->ndim, awhBiasParams->ndim);
 373     }
 374 }
 375
 376 /*! \brief Initializes the output file setup for the fricion output (note that the filename is not set here). */
 377 void OutputFile::initializeFrictionOutputFile(int                  subBlockStart,
 378                                               int                  numSubBlocks,
 379                                               const AwhBiasParams *awhBiasParams,
 380                                               EnergyUnit           energyUnit,
 381                                               real                 kTValue)
 382 {
 383     /* The first subblock with actual graph y-values is index 1 + ndim */
 384     numDim_               = awhBiasParams->ndim;
 385     int numTensorElements = (numDim_*(numDim_ + 1))/2;
 386
 387     /* The friction tensor elements are always the last subblocks */
 388     if (numSubBlocks < 1 + numDim_ + maxAwhGraphs + numTensorElements)
 389     {
 390         gmx_fatal(FARGS, "You requested friction tensor output, but the AWH data in the energy file does not contain the friction tensor");
 391     }
 392     GMX_ASSERT(numSubBlocks == 1 + numDim_ + maxAwhGraphs + numTensorElements, "The number of sub-blocks per bias should be 1 + ndim + maxAwhGraphs + (ndim*(ndim + 1))/2");
 393
 394     firstGraphSubBlock_ = subBlockStart + numSubBlocks - numTensorElements;
 395     numGraph_           = numTensorElements;
 396     useKTForEnergy_     = (energyUnit == EnergyUnit::KT);
 397     scaleFactor_.resize(numGraph_, useKTForEnergy_ ? 1 : kTValue);
 398     int numLegend       = numDim_ - 1 + numGraph_;
 399     legend_             = makeLegend(awhBiasParams, OutputFileType::Friction, numLegend);
 400     xLabel_             = "(nm or deg)";
 401     if (useKTForEnergy_)
 402     {
 403         yLabel_         = "friction/k\\sB\\NT (nm\\S-2\\Nps or rad\\S-2\\Nps)";
 404     }
 405     else
 406     {
 407         yLabel_         = "friction (kJ mol\\S-1\\Nnm\\S-2\\Nps or kJ mol\\S-1\\Nrad\\S-2\\Nps)";
 408     }
 409 }
 410
 411 AwhReader::AwhReader(const AwhParams  *awhParams,
 412                      int               numFileOptions,
 413                      const t_filenm   *filenames,
 414                      AwhGraphSelection awhGraphSelection,
 415                      EnergyUnit        energyUnit,
 416                      real              kT,
 417                      const t_enxblock *block) :
 418     kT_(kT)
 419 {
 420     GMX_RELEASE_ASSERT(block != nullptr, "NULL pointer passed to initAwhReader");
 421
 422     bool outputFriction = opt2bSet("-fric", numFileOptions, filenames);
 423
 424     /* The first subblock of each AWH block has metadata about
 425      * the number of subblocks belonging to that AWH block.
 426      * (block->nsub gives only the total number of subblocks and not how
 427      * they are distributed between the AWH biases).
 428      */
 429
 430     /* Keep track of the first subblock of this AWH */
 431     int subblockStart = 0;
 432     for (int k = 0; k < awhParams->numBias; k++)
 433     {
 434         AwhBiasParams              *awhBiasParams = &awhParams->awhBiasParams[k];
 435
 436         int                         numSubBlocks  = (int)block->sub[subblockStart].fval[0];
 437
 438         std::unique_ptr<OutputFile> awhOutputFile(new OutputFile(opt2fn("-o", numFileOptions, filenames), "AWH", awhParams->numBias, k));
 439
 440         awhOutputFile->initializeAwhOutputFile(subblockStart, numSubBlocks,
 441                                                awhBiasParams, awhGraphSelection,
 442                                                energyUnit, kT);
 443
 444         std::unique_ptr<OutputFile> frictionOutputFile;
 445         if (outputFriction)
 446         {
 447             frictionOutputFile = gmx::compat::make_unique<OutputFile>(opt2fn("-fric", numFileOptions, filenames), "Friction tensor", awhParams->numBias, k);
 448
 449             frictionOutputFile->initializeFrictionOutputFile(subblockStart, numSubBlocks, awhBiasParams, energyUnit, kT);
 450         }
 451
 452         biasOutputSetups_.emplace_back(BiasOutputSetup(subblockStart,
 453                                                        std::move(awhOutputFile),
 454                                                        std::move(frictionOutputFile)));
 455
 456         subblockStart += numSubBlocks;
 457     }
 458 }
 459
 460 FILE * OutputFile::openBiasOutputFile(double                  time,
 461                                       const gmx_output_env_t *oenv) const
 462 {
 463     std::string filename = baseFilename_ + gmx::formatString("_t%g.xvg", time);
 464
 465     FILE       *fp = xvgropen(filename.c_str(), title_.c_str(), xLabel_, yLabel_, oenv);
 466     xvgrLegend(fp, legend_, oenv);
 467
 468     return fp;
 469 }
 470
 471 /*! \brief Prints data selected by \p outputFile from \p block to \p fp */
 472 void OutputFile::writeData(const t_enxblock &block,
 473                            int               subBlockStart,
 474                            FILE             *fp) const
 475 {
 476     int numPoints = block.sub[subBlockStart + 1].nr;
 477     for (int j = 0; j < numPoints; j++)
 478     {
 479         /* Print the coordinates for numDim dimensions */
 480         for (int d = 0; d < numDim_; d++)
 481         {
 482             fprintf(fp, "  %8.4f", block.sub[subBlockStart + 1 + d].fval[j]);
 483         }
 484
 485         /* Print numGraph observables */
 486         for (int i = 0; i < numGraph_; i++)
 487         {
 488             fprintf(fp, "  %g", block.sub[firstGraphSubBlock_ + i].fval[j]*scaleFactor_[i]);
 489         }
 490
 491         fprintf(fp, "\n");
 492     }
 493 }
 494
 495 void AwhReader::processAwhFrame(const t_enxblock       &block,
 496                                 double                  time,
 497                                 const gmx_output_env_t *oenv) const
 498 {
 499     /* We look for AWH data every energy frame and count the no of AWH frames found. We only extract every 'skip' AWH frame. */
 500
 501     for (const auto &setup : biasOutputSetups_)
 502     {
 503         const int subStart = setup.subblockStart();
 504
 505         /* Each frame and AWH instance extracted generates one xvg file. */
 506         {
 507             const OutputFile &awhOutputFile = setup.awhOutputFile();
 508
 509             FILE             *fpAwh = awhOutputFile.openBiasOutputFile(time, oenv);
 510
 511             /* Now do the actual printing. Metadata in first subblock is treated separately. */
 512             fprintf(fpAwh, "# AWH metadata: target error = %.2f kT = %.2f kJ/mol\n",
 513                     block.sub[subStart].fval[1],
 514                     block.sub[subStart].fval[1]*kT_);
 515
 516             fprintf(fpAwh, "# AWH metadata: log sample weight = %4.2f\n",
 517                     block.sub[subStart].fval[2]);
 518
 519             awhOutputFile.writeData(block, subStart, fpAwh);
 520
 521             gmx_ffclose(fpAwh);
 522         }
 523
 524         const OutputFile *frictionOutputFile = setup.frictionOutputFile();
 525         if (frictionOutputFile != nullptr)
 526         {
 527             FILE *fpFriction = frictionOutputFile->openBiasOutputFile(time, oenv);
 528
 529             frictionOutputFile->writeData(block, subStart, fpFriction);
 530
 531             gmx_ffclose(fpFriction);
 532         }
 533     }
 534 }
 535
 536 /*! \brief The main function for the AWH tool */
 537 int gmx_awh(int argc, char *argv[])
 538 {
 539     const char         *desc[] = {
 540         "[THISMODULE] extracts AWH data from an energy file.",
 541         "One or two files are written per AWH bias per time frame.",
 542         "The bias index, if more than one, is appended to the file, as well as",
 543         "the time of the frame. By default only the PMF is printed.",
 544         "With [TT]-more[tt] the bias, target and coordinate distributions",
 545         "are also printed.",
 546         "With [TT]-more[tt] the bias, target and coordinate distributions",
 547         "are also printed, as well as the metric sqrt(det(friction_tensor))",
 548         "normalized such that the average is 1.",
 549         "Option [TT]-fric[tt] prints all components of the friction tensor",
 550         "to an additional set of files."
 551     };
 552     static gmx_bool     moreGraphs = FALSE;
 553     static int          skip       = 0;
 554     static gmx_bool     kTUnit     = FALSE;
 555     t_pargs             pa[]       = {
 556         { "-skip", FALSE, etINT,  {&skip},
 557           "Skip number of frames between data points" },
 558         { "-more", FALSE, etBOOL, {&moreGraphs},
 559           "Print more output" },
 560         { "-kt", FALSE, etBOOL, {&kTUnit},
 561           "Print free energy output in units of kT instead of kJ/mol" }
 562     };
 563
 564     ener_file_t         fp;
 565     t_inputrec          ir;
 566     gmx_enxnm_t        *enm = nullptr;
 567     t_enxframe         *frame;
 568     int                 nre;
 569     gmx_output_env_t   *oenv;
 570
 571     t_filenm            fnm[] = {
 572         { efEDR, "-f", nullptr,           ffREAD  },
 573         { efTPR, "-s", nullptr,           ffREAD  },
 574         { efXVG, "-o",    "awh",          ffWRITE },
 575         { efXVG, "-fric", "friction",     ffOPTWR }
 576     };
 577     const int           nfile  = asize(fnm);
 578     if (!parse_common_args(&argc, argv,
 579                            PCA_CAN_VIEW | PCA_CAN_BEGIN | PCA_CAN_END,
 580                            nfile, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
 581     {
 582         return 0;
 583     }
 584
 585     snew(frame, 1);
 586     fp = open_enx(ftp2fn(efEDR, nfile, fnm), "r");
 587     do_enxnms(fp, &nre, &enm);
 588
 589     /* We just need the AWH parameters from inputrec. These are used to initialize
 590        the AWH reader when we have a frame to read later on. */
 591     gmx_mtop_t mtop;
 592     int        natoms;
 593     matrix     box;
 594     read_tpx(ftp2fn(efTPR, nfile, fnm), &ir, box, &natoms, nullptr, nullptr, &mtop);
 595
 596     if (!ir.bDoAwh)
 597     {
 598         gmx_fatal(FARGS, "No AWH data in %s\n", opt2fn("-f", nfile, fnm));
 599     }
 600
 601     std::unique_ptr<AwhReader> awhReader;
 602
 603     /* Initiate counters */
 604     gmx_bool   haveFrame;
 605     int        awhFrameCounter = 0;
 606     int        timeCheck       = 0;
 607     do
 608     {
 609         haveFrame = do_enx(fp, frame);
 610
 611         bool        useFrame = false;
 612
 613         t_enxblock *block    = nullptr;
 614
 615         if (haveFrame)
 616         {
 617             timeCheck = check_times(frame->t);
 618
 619             if (timeCheck == 0)
 620             {
 621                 block = find_block_id_enxframe(frame, enxAWH, nullptr);
 622
 623                 if (block != nullptr)
 624                 {
 625                     if ((skip == 0) || (awhFrameCounter % skip == 0))
 626                     {
 627                         useFrame = true;
 628                     }
 629                     awhFrameCounter++;
 630                 }
 631             }
 632         }
 633
 634         if (useFrame)
 635         {
 636             /* We read a valid frame with an AWH block, so we can use it */
 637
 638             /* Currently we have the number of subblocks per AWH stored
 639              * in the first subblock (we cannot get this directly from the tpr),
 640              * so we have to read an AWH block before making the legends.
 641              */
 642             if (awhReader == nullptr)
 643             {
 644                 AwhGraphSelection awhGraphSelection = (moreGraphs ? AwhGraphSelection::All : AwhGraphSelection::Pmf);
 645                 EnergyUnit        energyUnit        = (kTUnit ? EnergyUnit::KT : EnergyUnit::KJPerMol);
 646                 awhReader =
 647                     std::unique_ptr<AwhReader>(new AwhReader(ir.awhParams,
 648                                                              nfile, fnm,
 649                                                              awhGraphSelection,
 650                                                              energyUnit, BOLTZ*ir.opts.ref_t[0],
 651                                                              block));
 652             }
 653
 654             awhReader->processAwhFrame(*block, frame->t, oenv);
 655         }
 656     }
 657     while (haveFrame && (timeCheck <= 0));
 658
 659     fprintf(stderr, "\n");
 660
 661     close_enx(fp);
 662
 663     return 0;
 664 }