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