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_;
 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 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_(false)
 322
 323 {
 324     baseFilename_ = filename.substr(0, filename.find('.'));
 325     title_        = baseTitle;
 326     if (numBias > 1)
 327     {
 328         baseFilename_ += gmx::formatString("%d", biasIndex + 1);
 329         title_        += gmx::formatString(" %d", biasIndex + 1);
 330     }
 331 }
 332
 333 void OutputFile::initializeAwhOutputFile(int                  subblockStart,
 334                                          int                  numSubBlocks,
 335                                          const AwhBiasParams *awhBiasParams,
 336                                          AwhGraphSelection    graphSelection,
 337                                          EnergyUnit           energyUnit,
 338                                          real                 kTValue)
 339 {
 340     /* The first subblock with actual graph y-values is index 1 + ndim */
 341     numDim_             = awhBiasParams->ndim;
 342     firstGraphSubBlock_ = subblockStart + 1 + numDim_;
 343     if (graphSelection == AwhGraphSelection::All)
 344     {
 345         /* There are one metadata and ndim coordinate blocks */
 346         numGraph_       = std::min(numSubBlocks - 1 - numDim_,
 347                                    maxAwhGraphs);
 348     }
 349     else
 350     {
 351         numGraph_       = 1;
 352     }
 353     useKTForEnergy_     = (energyUnit == EnergyUnit::KT);
 354     scaleFactor_.resize(numGraph_, 1);
 355     if (!useKTForEnergy_)
 356     {
 357         /* The first two graphs are in units of energy, multiply by kT */
 358         std::fill(scaleFactor_.begin(), scaleFactor_.begin() + std::min(2, numGraph_), kTValue);
 359     }
 360     int numLegend = numDim_ - 1 + numGraph_;
 361     legend_       = makeLegend(awhBiasParams, OutputFileType::Awh, numLegend);
 362     /* We could have both length and angle coordinates in a single bias */
 363     xLabel_       = "(nm or deg)";
 364     yLabel_       = useKTForEnergy_ ? "(k\\sB\\NT)" : "(kJ/mol)";
 365     if (graphSelection == AwhGraphSelection::All)
 366     {
 367         yLabel_ += gmx::formatString(", (nm\\S-%d\\N or rad\\S-%d\\N), (-)",
 368                                      awhBiasParams->ndim, awhBiasParams->ndim);
 369     }
 370 }
 371
 372 /*! \brief Initializes the output file setup for the fricion output (note that the filename is not set here). */
 373 void OutputFile::initializeFrictionOutputFile(int                  subBlockStart,
 374                                               int                  numSubBlocks,
 375                                               const AwhBiasParams *awhBiasParams,
 376                                               EnergyUnit           energyUnit,
 377                                               real                 kTValue)
 378 {
 379     /* The first subblock with actual graph y-values is index 1 + ndim */
 380     numDim_               = awhBiasParams->ndim;
 381     int numTensorElements = (numDim_*(numDim_ + 1))/2;
 382
 383     /* The friction tensor elements are always the last subblocks */
 384     if (numSubBlocks < 1 + numDim_ + maxAwhGraphs + numTensorElements)
 385     {
 386         gmx_fatal(FARGS, "You requested friction tensor output, but the AWH data in the energy file does not contain the friction tensor");
 387     }
 388     GMX_ASSERT(numSubBlocks == 1 + numDim_ + maxAwhGraphs + numTensorElements, "The number of sub-blocks per bias should be 1 + ndim + maxAwhGraphs + (ndim*(ndim + 1))/2");
 389
 390     firstGraphSubBlock_ = subBlockStart + numSubBlocks - numTensorElements;
 391     numGraph_           = numTensorElements;
 392     useKTForEnergy_     = (energyUnit == EnergyUnit::KT);
 393     scaleFactor_.resize(numGraph_, useKTForEnergy_ ? 1 : kTValue);
 394     int numLegend       = numDim_ - 1 + numGraph_;
 395     legend_             = makeLegend(awhBiasParams, OutputFileType::Friction, numLegend);
 396     xLabel_             = "(nm or deg)";
 397     if (useKTForEnergy_)
 398     {
 399         yLabel_         = "friction/k\\sB\\NT (nm\\S-2\\Nps or rad\\S-2\\Nps)";
 400     }
 401     else
 402     {
 403         yLabel_         = "friction (kJ mol\\S-1\\Nnm\\S-2\\Nps or kJ mol\\S-1\\Nrad\\S-2\\Nps)";
 404     }
 405 }
 406
 407 AwhReader::AwhReader(const AwhParams  *awhParams,
 408                      int               numFileOptions,
 409                      const t_filenm   *filenames,
 410                      AwhGraphSelection awhGraphSelection,
 411                      EnergyUnit        energyUnit,
 412                      real              kT,
 413                      const t_enxblock *block) :
 414     kT_(kT)
 415 {
 416     GMX_RELEASE_ASSERT(block != nullptr, "NULL pointer passed to initAwhReader");
 417
 418     bool outputFriction = opt2bSet("-fric", numFileOptions, filenames);
 419
 420     /* The first subblock of each AWH block has metadata about
 421      * the number of subblocks belonging to that AWH block.
 422      * (block->nsub gives only the total number of subblocks and not how
 423      * they are distributed between the AWH biases).
 424      */
 425
 426     /* Keep track of the first subblock of this AWH */
 427     int subblockStart = 0;
 428     for (int k = 0; k < awhParams->numBias; k++)
 429     {
 430         AwhBiasParams              *awhBiasParams = &awhParams->awhBiasParams[k];
 431
 432         int                         numSubBlocks  = static_cast<int>(block->sub[subblockStart].fval[0]);
 433
 434         std::unique_ptr<OutputFile> awhOutputFile(new OutputFile(opt2fn("-o", numFileOptions, filenames), "AWH", awhParams->numBias, k));
 435
 436         awhOutputFile->initializeAwhOutputFile(subblockStart, numSubBlocks,
 437                                                awhBiasParams, awhGraphSelection,
 438                                                energyUnit, kT);
 439
 440         std::unique_ptr<OutputFile> frictionOutputFile;
 441         if (outputFriction)
 442         {
 443             frictionOutputFile = gmx::compat::make_unique<OutputFile>(opt2fn("-fric", numFileOptions, filenames), "Friction tensor", awhParams->numBias, k);
 444
 445             frictionOutputFile->initializeFrictionOutputFile(subblockStart, numSubBlocks, awhBiasParams, energyUnit, kT);
 446         }
 447
 448         biasOutputSetups_.emplace_back(BiasOutputSetup(subblockStart,
 449                                                        std::move(awhOutputFile),
 450                                                        std::move(frictionOutputFile)));
 451
 452         subblockStart += numSubBlocks;
 453     }
 454 }
 455
 456 FILE * OutputFile::openBiasOutputFile(double                  time,
 457                                       const gmx_output_env_t *oenv) const
 458 {
 459     std::string filename = baseFilename_ + gmx::formatString("_t%g.xvg", time);
 460
 461     FILE       *fp = xvgropen(filename.c_str(), title_.c_str(), xLabel_, yLabel_, oenv);
 462     xvgrLegend(fp, legend_, oenv);
 463
 464     return fp;
 465 }
 466
 467 /*! \brief Prints data selected by \p outputFile from \p block to \p fp */
 468 void OutputFile::writeData(const t_enxblock &block,
 469                            int               subBlockStart,
 470                            FILE             *fp) const
 471 {
 472     int numPoints = block.sub[subBlockStart + 1].nr;
 473     for (int j = 0; j < numPoints; j++)
 474     {
 475         /* Print the coordinates for numDim dimensions */
 476         for (int d = 0; d < numDim_; d++)
 477         {
 478             fprintf(fp, "  %8.4f", block.sub[subBlockStart + 1 + d].fval[j]);
 479         }
 480
 481         /* Print numGraph observables */
 482         for (int i = 0; i < numGraph_; i++)
 483         {
 484             fprintf(fp, "  %g", block.sub[firstGraphSubBlock_ + i].fval[j]*scaleFactor_[i]);
 485         }
 486
 487         fprintf(fp, "\n");
 488     }
 489 }
 490
 491 void AwhReader::processAwhFrame(const t_enxblock       &block,
 492                                 double                  time,
 493                                 const gmx_output_env_t *oenv) const
 494 {
 495     /* We look for AWH data every energy frame and count the no of AWH frames found. We only extract every 'skip' AWH frame. */
 496
 497     for (const auto &setup : biasOutputSetups_)
 498     {
 499         const int subStart = setup.subblockStart();
 500
 501         /* Each frame and AWH instance extracted generates one xvg file. */
 502         {
 503             const OutputFile &awhOutputFile = setup.awhOutputFile();
 504
 505             FILE             *fpAwh = awhOutputFile.openBiasOutputFile(time, oenv);
 506
 507             /* Now do the actual printing. Metadata in first subblock is treated separately. */
 508             fprintf(fpAwh, "# AWH metadata: target error = %.2f kT = %.2f kJ/mol\n",
 509                     block.sub[subStart].fval[1],
 510                     block.sub[subStart].fval[1]*kT_);
 511
 512             fprintf(fpAwh, "# AWH metadata: log sample weight = %4.2f\n",
 513                     block.sub[subStart].fval[2]);
 514
 515             awhOutputFile.writeData(block, subStart, fpAwh);
 516
 517             gmx_ffclose(fpAwh);
 518         }
 519
 520         const OutputFile *frictionOutputFile = setup.frictionOutputFile();
 521         if (frictionOutputFile != nullptr)
 522         {
 523             FILE *fpFriction = frictionOutputFile->openBiasOutputFile(time, oenv);
 524
 525             frictionOutputFile->writeData(block, subStart, fpFriction);
 526
 527             gmx_ffclose(fpFriction);
 528         }
 529     }
 530 }
 531
 532 /*! \brief The main function for the AWH tool */
 533 int gmx_awh(int argc, char *argv[])
 534 {
 535     const char         *desc[] = {
 536         "[THISMODULE] extracts AWH data from an energy file.",
 537         "One or two files are written per AWH bias per time frame.",
 538         "The bias index, if more than one, is appended to the file, as well as",
 539         "the time of the frame. By default only the PMF is printed.",
 540         "With [TT]-more[tt] the bias, target and coordinate distributions",
 541         "are also printed.",
 542         "With [TT]-more[tt] the bias, target and coordinate distributions",
 543         "are also printed, as well as the metric sqrt(det(friction_tensor))",
 544         "normalized such that the average is 1.",
 545         "Option [TT]-fric[tt] prints all components of the friction tensor",
 546         "to an additional set of files."
 547     };
 548     static gmx_bool     moreGraphs = FALSE;
 549     static int          skip       = 0;
 550     static gmx_bool     kTUnit     = FALSE;
 551     t_pargs             pa[]       = {
 552         { "-skip", FALSE, etINT,  {&skip},
 553           "Skip number of frames between data points" },
 554         { "-more", FALSE, etBOOL, {&moreGraphs},
 555           "Print more output" },
 556         { "-kt", FALSE, etBOOL, {&kTUnit},
 557           "Print free energy output in units of kT instead of kJ/mol" }
 558     };
 559
 560     ener_file_t         fp;
 561     t_inputrec          ir;
 562     gmx_enxnm_t        *enm = nullptr;
 563     t_enxframe         *frame;
 564     int                 nre;
 565     gmx_output_env_t   *oenv;
 566
 567     t_filenm            fnm[] = {
 568         { efEDR, "-f", nullptr,           ffREAD  },
 569         { efTPR, "-s", nullptr,           ffREAD  },
 570         { efXVG, "-o",    "awh",          ffWRITE },
 571         { efXVG, "-fric", "friction",     ffOPTWR }
 572     };
 573     const int           nfile  = asize(fnm);
 574     if (!parse_common_args(&argc, argv,
 575                            PCA_CAN_VIEW | PCA_CAN_BEGIN | PCA_CAN_END,
 576                            nfile, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
 577     {
 578         return 0;
 579     }
 580
 581     snew(frame, 1);
 582     fp = open_enx(ftp2fn(efEDR, nfile, fnm), "r");
 583     do_enxnms(fp, &nre, &enm);
 584
 585     /* We just need the AWH parameters from inputrec. These are used to initialize
 586        the AWH reader when we have a frame to read later on. */
 587     gmx_mtop_t mtop;
 588     int        natoms;
 589     matrix     box;
 590     read_tpx(ftp2fn(efTPR, nfile, fnm), &ir, box, &natoms, nullptr, nullptr, &mtop);
 591
 592     if (!ir.bDoAwh)
 593     {
 594         gmx_fatal(FARGS, "No AWH data in %s\n", opt2fn("-f", nfile, fnm));
 595     }
 596
 597     std::unique_ptr<AwhReader> awhReader;
 598
 599     /* Initiate counters */
 600     gmx_bool   haveFrame;
 601     int        awhFrameCounter = 0;
 602     int        timeCheck       = 0;
 603     do
 604     {
 605         haveFrame = do_enx(fp, frame);
 606
 607         bool        useFrame = false;
 608
 609         t_enxblock *block    = nullptr;
 610
 611         if (haveFrame)
 612         {
 613             timeCheck = check_times(frame->t);
 614
 615             if (timeCheck == 0)
 616             {
 617                 block = find_block_id_enxframe(frame, enxAWH, nullptr);
 618
 619                 if (block != nullptr)
 620                 {
 621                     if ((skip == 0) || (awhFrameCounter % skip == 0))
 622                     {
 623                         useFrame = true;
 624                     }
 625                     awhFrameCounter++;
 626                 }
 627             }
 628         }
 629
 630         if (useFrame)
 631         {
 632             /* We read a valid frame with an AWH block, so we can use it */
 633
 634             /* Currently we have the number of subblocks per AWH stored
 635              * in the first subblock (we cannot get this directly from the tpr),
 636              * so we have to read an AWH block before making the legends.
 637              */
 638             if (awhReader == nullptr)
 639             {
 640                 AwhGraphSelection awhGraphSelection = (moreGraphs ? AwhGraphSelection::All : AwhGraphSelection::Pmf);
 641                 EnergyUnit        energyUnit        = (kTUnit ? EnergyUnit::KT : EnergyUnit::KJPerMol);
 642                 awhReader =
 643                     gmx::compat::make_unique<AwhReader>(ir.awhParams,
 644                                                         nfile, fnm,
 645                                                         awhGraphSelection,
 646                                                         energyUnit, BOLTZ*ir.opts.ref_t[0],
 647                                                         block);
 648             }
 649
 650             awhReader->processAwhFrame(*block, frame->t, oenv);
 651         }
 652     }
 653     while (haveFrame && (timeCheck <= 0));
 654
 655     fprintf(stderr, "\n");
 656
 657     close_enx(fp);
 658
 659     return 0;
 660 }