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