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