2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
38 * Tool for extracting AWH (Accelerated Weight Histogram) data from energy files.
40 * \author Viveca Lindahl
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"
78 using gmx::AwhBiasParams;
83 //! Enum for choosing AWH graph output
84 enum class AwhGraphSelection
86 Pmf, //!< Only the PMF
87 All //!< All possible AWH graphs
90 //! Energy unit options
99 /*! \brief All meta-data that is shared for one output file type for one bias */
103 /*! \brief Constructor, Set the output base file name and title.
105 * Result is a valid object, but will produce empty output files.
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.
112 OutputFile(const std::string &filename,
113 const std::string &baseTitle,
117 /*! \brief Initializes the output file setup for the AWH output.
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.
126 void initializeAwhOutputFile(int subBlockStart,
128 const AwhBiasParams *awhBiasParams,
129 AwhGraphSelection graphSelection,
130 EnergyUnit energyUnit,
133 /*! \brief Initializes the output file setup for the fricion output.
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.
141 void initializeFrictionOutputFile(int subBlockStart,
143 const AwhBiasParams *awhBiasParams,
144 EnergyUnit energyUnit,
147 /*! \brief Opens a single output file for a bias, prints title and legends.
149 * \param[in] time The time for of frame to be written.
150 * \param[in] oenv The output environment.
151 * \returns the file pointer.
153 FILE * openBiasOutputFile(double time,
154 const gmx_output_env_t *oenv) const;
156 /*! \brief Writes data selected from \p block to file.
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.
162 void writeData(const t_enxblock &block,
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. */
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. */
180 /*! \brief All meta-data that is shared for all output files for one bias */
181 struct BiasOutputSetup
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))
193 //! Return the AWH output file data.
194 const OutputFile &awhOutputFile() const
196 GMX_RELEASE_ASSERT(awhOutputFile_ != nullptr, "awhOutputFile() called without initialized AWH output file");
198 return *awhOutputFile_;
201 //! Return the a pointer to the friction output file data, can return nullptr
202 const OutputFile *frictionOutputFile() const
204 return frictionOutputFile_.get();
207 //! Return the starting subblock.
208 int subblockStart() const
210 return subblockStart_;
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 */
219 /*! \brief All options and meta-data needed for the AWH output */
224 AwhReader(const AwhParams *awhParams,
226 const t_filenm *filenames,
227 AwhGraphSelection awhGraphSelection,
228 EnergyUnit energyUnit,
230 const t_enxblock *block);
232 //! Extract and print AWH data for an AWH block of one energy frame
233 void processAwhFrame(const t_enxblock &block,
235 const gmx_output_env_t *oenv) const;
238 std::vector<BiasOutputSetup> biasOutputSetups_; /**< The output setups, one for each AWH bias. */
240 const real kT_; /**< kB*T in kJ/mol. */
246 //! Enum for selecting output file type.
247 enum class OutputFileType
249 Awh, //!< AWH, all data except friction tensor.
250 Friction //!< The full friction tensor.
253 /*! \brief The maximum number of observables per subblock, not including the full friction tensor.
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.
259 constexpr int maxAwhGraphs = 6;
261 /*! \brief Constructs a legend for a standard awh output file */
262 std::vector<std::string>makeLegend(const AwhBiasParams *awhBiasParams,
263 OutputFileType outputFileType,
266 const std::array<std::string, maxAwhGraphs> legendBase =
272 "Target ref value distr",
276 std::vector<std::string> legend;
277 /* Give legends to dimensions higher than the first */
278 for (int d = 1; d < awhBiasParams->ndim; d++)
280 legend.push_back(gmx::formatString("dim%d", d));
283 switch (outputFileType)
285 case OutputFileType::Awh:
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())
291 legend.push_back(legendBase[legendBaseIndex]);
296 case OutputFileType::Friction:
297 for (int i0 = 0; i0 < awhBiasParams->ndim; i0++)
299 for (int i1 = 0; i1 <= i0; i1++)
301 legend.push_back(gmx::formatString("%d,%d", i0, i1));
307 GMX_RELEASE_ASSERT(legend.size() == numLegend, "The number of legends requested for printing and the number generated should be the same");
314 OutputFile::OutputFile(const std::string &filename,
315 const std::string &baseTitle,
319 firstGraphSubBlock_(0),
324 // cppcheck-suppress useInitializationList
325 baseFilename_ = filename.substr(0, filename.find('.'));
329 baseFilename_ += gmx::formatString("%d", biasIndex + 1);
330 title_ += gmx::formatString(" %d", biasIndex + 1);
334 void OutputFile::initializeAwhOutputFile(int subblockStart,
336 const AwhBiasParams *awhBiasParams,
337 AwhGraphSelection graphSelection,
338 EnergyUnit energyUnit,
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)
346 /* There are one metadata and ndim coordinate blocks */
347 numGraph_ = std::min(numSubBlocks - 1 - numDim_,
354 useKTForEnergy_ = (energyUnit == EnergyUnit::KT);
355 scaleFactor_.resize(numGraph_, 1);
356 if (!useKTForEnergy_)
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);
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)
368 yLabel_ += gmx::formatString(", (nm\\S-%d\\N or rad\\S-%d\\N), (-)",
369 awhBiasParams->ndim, awhBiasParams->ndim);
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,
376 const AwhBiasParams *awhBiasParams,
377 EnergyUnit energyUnit,
380 /* The first subblock with actual graph y-values is index 1 + ndim */
381 numDim_ = awhBiasParams->ndim;
382 int numTensorElements = (numDim_*(numDim_ + 1))/2;
384 /* The friction tensor elements are always the last subblocks */
385 if (numSubBlocks < 1 + numDim_ + maxAwhGraphs + numTensorElements)
387 gmx_fatal(FARGS, "You requested friction tensor output, but the AWH data in the energy file does not contain the friction tensor");
389 GMX_ASSERT(numSubBlocks == 1 + numDim_ + maxAwhGraphs + numTensorElements, "The number of sub-blocks per bias should be 1 + ndim + maxAwhGraphs + (ndim*(ndim + 1))/2");
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)";
400 yLabel_ = "friction/k\\sB\\NT (nm\\S-2\\Nps or rad\\S-2\\Nps)";
404 yLabel_ = "friction (kJ mol\\S-1\\Nnm\\S-2\\Nps or kJ mol\\S-1\\Nrad\\S-2\\Nps)";
408 AwhReader::AwhReader(const AwhParams *awhParams,
410 const t_filenm *filenames,
411 AwhGraphSelection awhGraphSelection,
412 EnergyUnit energyUnit,
414 const t_enxblock *block) :
417 GMX_RELEASE_ASSERT(block != nullptr, "NULL pointer passed to initAwhReader");
419 bool outputFriction = opt2bSet("-fric", numFileOptions, filenames);
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).
427 /* Keep track of the first subblock of this AWH */
428 int subblockStart = 0;
429 for (int k = 0; k < awhParams->numBias; k++)
431 AwhBiasParams *awhBiasParams = &awhParams->awhBiasParams[k];
433 int numSubBlocks = (int)block->sub[subblockStart].fval[0];
435 std::unique_ptr<OutputFile> awhOutputFile(new OutputFile(opt2fn("-o", numFileOptions, filenames), "AWH", awhParams->numBias, k));
437 awhOutputFile->initializeAwhOutputFile(subblockStart, numSubBlocks,
438 awhBiasParams, awhGraphSelection,
441 std::unique_ptr<OutputFile> frictionOutputFile;
444 frictionOutputFile = gmx::compat::make_unique<OutputFile>(opt2fn("-fric", numFileOptions, filenames), "Friction tensor", awhParams->numBias, k);
446 frictionOutputFile->initializeFrictionOutputFile(subblockStart, numSubBlocks, awhBiasParams, energyUnit, kT);
449 biasOutputSetups_.emplace_back(BiasOutputSetup(subblockStart,
450 std::move(awhOutputFile),
451 std::move(frictionOutputFile)));
453 subblockStart += numSubBlocks;
457 FILE * OutputFile::openBiasOutputFile(double time,
458 const gmx_output_env_t *oenv) const
460 std::string filename = baseFilename_ + gmx::formatString("_t%g.xvg", time);
462 FILE *fp = xvgropen(filename.c_str(), title_.c_str(), xLabel_, yLabel_, oenv);
463 xvgrLegend(fp, legend_, oenv);
468 /*! \brief Prints data selected by \p outputFile from \p block to \p fp */
469 void OutputFile::writeData(const t_enxblock &block,
473 int numPoints = block.sub[subBlockStart + 1].nr;
474 for (int j = 0; j < numPoints; j++)
476 /* Print the coordinates for numDim dimensions */
477 for (int d = 0; d < numDim_; d++)
479 fprintf(fp, " %8.4f", block.sub[subBlockStart + 1 + d].fval[j]);
482 /* Print numGraph observables */
483 for (int i = 0; i < numGraph_; i++)
485 fprintf(fp, " %g", block.sub[firstGraphSubBlock_ + i].fval[j]*scaleFactor_[i]);
492 void AwhReader::processAwhFrame(const t_enxblock &block,
494 const gmx_output_env_t *oenv) const
496 /* We look for AWH data every energy frame and count the no of AWH frames found. We only extract every 'skip' AWH frame. */
498 for (const auto &setup : biasOutputSetups_)
500 const int subStart = setup.subblockStart();
502 /* Each frame and AWH instance extracted generates one xvg file. */
504 const OutputFile &awhOutputFile = setup.awhOutputFile();
506 FILE *fpAwh = awhOutputFile.openBiasOutputFile(time, oenv);
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_);
513 fprintf(fpAwh, "# AWH metadata: log sample weight = %4.2f\n",
514 block.sub[subStart].fval[2]);
516 awhOutputFile.writeData(block, subStart, fpAwh);
521 const OutputFile *frictionOutputFile = setup.frictionOutputFile();
522 if (frictionOutputFile != nullptr)
524 FILE *fpFriction = frictionOutputFile->openBiasOutputFile(time, oenv);
526 frictionOutputFile->writeData(block, subStart, fpFriction);
528 gmx_ffclose(fpFriction);
533 /*! \brief The main function for the AWH tool */
534 int gmx_awh(int argc, char *argv[])
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",
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."
549 static gmx_bool moreGraphs = FALSE;
551 static gmx_bool kTUnit = FALSE;
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" }
563 gmx_enxnm_t *enm = nullptr;
566 gmx_output_env_t *oenv;
569 { efEDR, "-f", nullptr, ffREAD },
570 { efTPR, "-s", nullptr, ffREAD },
571 { efXVG, "-o", "awh", ffWRITE },
572 { efXVG, "-fric", "friction", ffOPTWR }
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))
583 fp = open_enx(ftp2fn(efEDR, nfile, fnm), "r");
584 do_enxnms(fp, &nre, &enm);
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. */
591 read_tpx(ftp2fn(efTPR, nfile, fnm), &ir, box, &natoms, nullptr, nullptr, &mtop);
595 gmx_fatal(FARGS, "No AWH data in %s\n", opt2fn("-f", nfile, fnm));
598 std::unique_ptr<AwhReader> awhReader;
600 /* Initiate counters */
602 int awhFrameCounter = 0;
606 haveFrame = do_enx(fp, frame);
608 bool useFrame = false;
610 t_enxblock *block = nullptr;
614 timeCheck = check_times(frame->t);
618 block = find_block_id_enxframe(frame, enxAWH, nullptr);
620 if (block != nullptr)
622 if ((skip == 0) || (awhFrameCounter % skip == 0))
633 /* We read a valid frame with an AWH block, so we can use it */
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.
639 if (awhReader == nullptr)
641 AwhGraphSelection awhGraphSelection = (moreGraphs ? AwhGraphSelection::All : AwhGraphSelection::Pmf);
642 EnergyUnit energyUnit = (kTUnit ? EnergyUnit::KT : EnergyUnit::KJPerMol);
644 gmx::compat::make_unique<AwhReader>(ir.awhParams,
647 energyUnit, BOLTZ*ir.opts.ref_t[0],
651 awhReader->processAwhFrame(*block, frame->t, oenv);
654 while (haveFrame && (timeCheck <= 0));
656 fprintf(stderr, "\n");