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