src/gromacs/gmxana/gmx_awh.cpp

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