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