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),
321 useKTForEnergy_(false)
324 baseFilename_ = filename.substr(0, filename.find('.'));
328 baseFilename_ += gmx::formatString("%d", biasIndex + 1);
329 title_ += gmx::formatString(" %d", biasIndex + 1);
333 void OutputFile::initializeAwhOutputFile(int subblockStart,
335 const AwhBiasParams *awhBiasParams,
336 AwhGraphSelection graphSelection,
337 EnergyUnit energyUnit,
340 /* The first subblock with actual graph y-values is index 1 + ndim */
341 numDim_ = awhBiasParams->ndim;
342 firstGraphSubBlock_ = subblockStart + 1 + numDim_;
343 if (graphSelection == AwhGraphSelection::All)
345 /* There are one metadata and ndim coordinate blocks */
346 numGraph_ = std::min(numSubBlocks - 1 - numDim_,
353 useKTForEnergy_ = (energyUnit == EnergyUnit::KT);
354 scaleFactor_.resize(numGraph_, 1);
355 if (!useKTForEnergy_)
357 /* The first two graphs are in units of energy, multiply by kT */
358 std::fill(scaleFactor_.begin(), scaleFactor_.begin() + std::min(2, numGraph_), kTValue);
360 int numLegend = numDim_ - 1 + numGraph_;
361 legend_ = makeLegend(awhBiasParams, OutputFileType::Awh, numLegend);
362 /* We could have both length and angle coordinates in a single bias */
363 xLabel_ = "(nm or deg)";
364 yLabel_ = useKTForEnergy_ ? "(k\\sB\\NT)" : "(kJ/mol)";
365 if (graphSelection == AwhGraphSelection::All)
367 yLabel_ += gmx::formatString(", (nm\\S-%d\\N or rad\\S-%d\\N), (-)",
368 awhBiasParams->ndim, awhBiasParams->ndim);
372 /*! \brief Initializes the output file setup for the fricion output (note that the filename is not set here). */
373 void OutputFile::initializeFrictionOutputFile(int subBlockStart,
375 const AwhBiasParams *awhBiasParams,
376 EnergyUnit energyUnit,
379 /* The first subblock with actual graph y-values is index 1 + ndim */
380 numDim_ = awhBiasParams->ndim;
381 int numTensorElements = (numDim_*(numDim_ + 1))/2;
383 /* The friction tensor elements are always the last subblocks */
384 if (numSubBlocks < 1 + numDim_ + maxAwhGraphs + numTensorElements)
386 gmx_fatal(FARGS, "You requested friction tensor output, but the AWH data in the energy file does not contain the friction tensor");
388 GMX_ASSERT(numSubBlocks == 1 + numDim_ + maxAwhGraphs + numTensorElements, "The number of sub-blocks per bias should be 1 + ndim + maxAwhGraphs + (ndim*(ndim + 1))/2");
390 firstGraphSubBlock_ = subBlockStart + numSubBlocks - numTensorElements;
391 numGraph_ = numTensorElements;
392 useKTForEnergy_ = (energyUnit == EnergyUnit::KT);
393 scaleFactor_.resize(numGraph_, useKTForEnergy_ ? 1 : kTValue);
394 int numLegend = numDim_ - 1 + numGraph_;
395 legend_ = makeLegend(awhBiasParams, OutputFileType::Friction, numLegend);
396 xLabel_ = "(nm or deg)";
399 yLabel_ = "friction/k\\sB\\NT (nm\\S-2\\Nps or rad\\S-2\\Nps)";
403 yLabel_ = "friction (kJ mol\\S-1\\Nnm\\S-2\\Nps or kJ mol\\S-1\\Nrad\\S-2\\Nps)";
407 AwhReader::AwhReader(const AwhParams *awhParams,
409 const t_filenm *filenames,
410 AwhGraphSelection awhGraphSelection,
411 EnergyUnit energyUnit,
413 const t_enxblock *block) :
416 GMX_RELEASE_ASSERT(block != nullptr, "NULL pointer passed to initAwhReader");
418 bool outputFriction = opt2bSet("-fric", numFileOptions, filenames);
420 /* The first subblock of each AWH block has metadata about
421 * the number of subblocks belonging to that AWH block.
422 * (block->nsub gives only the total number of subblocks and not how
423 * they are distributed between the AWH biases).
426 /* Keep track of the first subblock of this AWH */
427 int subblockStart = 0;
428 for (int k = 0; k < awhParams->numBias; k++)
430 AwhBiasParams *awhBiasParams = &awhParams->awhBiasParams[k];
432 int numSubBlocks = static_cast<int>(block->sub[subblockStart].fval[0]);
434 std::unique_ptr<OutputFile> awhOutputFile(new OutputFile(opt2fn("-o", numFileOptions, filenames), "AWH", awhParams->numBias, k));
436 awhOutputFile->initializeAwhOutputFile(subblockStart, numSubBlocks,
437 awhBiasParams, awhGraphSelection,
440 std::unique_ptr<OutputFile> frictionOutputFile;
443 frictionOutputFile = gmx::compat::make_unique<OutputFile>(opt2fn("-fric", numFileOptions, filenames), "Friction tensor", awhParams->numBias, k);
445 frictionOutputFile->initializeFrictionOutputFile(subblockStart, numSubBlocks, awhBiasParams, energyUnit, kT);
448 biasOutputSetups_.emplace_back(BiasOutputSetup(subblockStart,
449 std::move(awhOutputFile),
450 std::move(frictionOutputFile)));
452 subblockStart += numSubBlocks;
456 FILE * OutputFile::openBiasOutputFile(double time,
457 const gmx_output_env_t *oenv) const
459 std::string filename = baseFilename_ + gmx::formatString("_t%g.xvg", time);
461 FILE *fp = xvgropen(filename.c_str(), title_.c_str(), xLabel_, yLabel_, oenv);
462 xvgrLegend(fp, legend_, oenv);
467 /*! \brief Prints data selected by \p outputFile from \p block to \p fp */
468 void OutputFile::writeData(const t_enxblock &block,
472 int numPoints = block.sub[subBlockStart + 1].nr;
473 for (int j = 0; j < numPoints; j++)
475 /* Print the coordinates for numDim dimensions */
476 for (int d = 0; d < numDim_; d++)
478 fprintf(fp, " %8.4f", block.sub[subBlockStart + 1 + d].fval[j]);
481 /* Print numGraph observables */
482 for (int i = 0; i < numGraph_; i++)
484 fprintf(fp, " %g", block.sub[firstGraphSubBlock_ + i].fval[j]*scaleFactor_[i]);
491 void AwhReader::processAwhFrame(const t_enxblock &block,
493 const gmx_output_env_t *oenv) const
495 /* We look for AWH data every energy frame and count the no of AWH frames found. We only extract every 'skip' AWH frame. */
497 for (const auto &setup : biasOutputSetups_)
499 const int subStart = setup.subblockStart();
501 /* Each frame and AWH instance extracted generates one xvg file. */
503 const OutputFile &awhOutputFile = setup.awhOutputFile();
505 FILE *fpAwh = awhOutputFile.openBiasOutputFile(time, oenv);
507 /* Now do the actual printing. Metadata in first subblock is treated separately. */
508 fprintf(fpAwh, "# AWH metadata: target error = %.2f kT = %.2f kJ/mol\n",
509 block.sub[subStart].fval[1],
510 block.sub[subStart].fval[1]*kT_);
512 fprintf(fpAwh, "# AWH metadata: log sample weight = %4.2f\n",
513 block.sub[subStart].fval[2]);
515 awhOutputFile.writeData(block, subStart, fpAwh);
520 const OutputFile *frictionOutputFile = setup.frictionOutputFile();
521 if (frictionOutputFile != nullptr)
523 FILE *fpFriction = frictionOutputFile->openBiasOutputFile(time, oenv);
525 frictionOutputFile->writeData(block, subStart, fpFriction);
527 gmx_ffclose(fpFriction);
532 /*! \brief The main function for the AWH tool */
533 int gmx_awh(int argc, char *argv[])
535 const char *desc[] = {
536 "[THISMODULE] extracts AWH data from an energy file.",
537 "One or two files are written per AWH bias per time frame.",
538 "The bias index, if more than one, is appended to the file, as well as",
539 "the time of the frame. By default only the PMF is printed.",
540 "With [TT]-more[tt] the bias, target and coordinate distributions",
542 "With [TT]-more[tt] the bias, target and coordinate distributions",
543 "are also printed, as well as the metric sqrt(det(friction_tensor))",
544 "normalized such that the average is 1.",
545 "Option [TT]-fric[tt] prints all components of the friction tensor",
546 "to an additional set of files."
548 static gmx_bool moreGraphs = FALSE;
550 static gmx_bool kTUnit = FALSE;
552 { "-skip", FALSE, etINT, {&skip},
553 "Skip number of frames between data points" },
554 { "-more", FALSE, etBOOL, {&moreGraphs},
555 "Print more output" },
556 { "-kt", FALSE, etBOOL, {&kTUnit},
557 "Print free energy output in units of kT instead of kJ/mol" }
562 gmx_enxnm_t *enm = nullptr;
565 gmx_output_env_t *oenv;
568 { efEDR, "-f", nullptr, ffREAD },
569 { efTPR, "-s", nullptr, ffREAD },
570 { efXVG, "-o", "awh", ffWRITE },
571 { efXVG, "-fric", "friction", ffOPTWR }
573 const int nfile = asize(fnm);
574 if (!parse_common_args(&argc, argv,
575 PCA_CAN_VIEW | PCA_CAN_BEGIN | PCA_CAN_END,
576 nfile, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
582 fp = open_enx(ftp2fn(efEDR, nfile, fnm), "r");
583 do_enxnms(fp, &nre, &enm);
585 /* We just need the AWH parameters from inputrec. These are used to initialize
586 the AWH reader when we have a frame to read later on. */
590 read_tpx(ftp2fn(efTPR, nfile, fnm), &ir, box, &natoms, nullptr, nullptr, &mtop);
594 gmx_fatal(FARGS, "No AWH data in %s\n", opt2fn("-f", nfile, fnm));
597 std::unique_ptr<AwhReader> awhReader;
599 /* Initiate counters */
601 int awhFrameCounter = 0;
605 haveFrame = do_enx(fp, frame);
607 bool useFrame = false;
609 t_enxblock *block = nullptr;
613 timeCheck = check_times(frame->t);
617 block = find_block_id_enxframe(frame, enxAWH, nullptr);
619 if (block != nullptr)
621 if ((skip == 0) || (awhFrameCounter % skip == 0))
632 /* We read a valid frame with an AWH block, so we can use it */
634 /* Currently we have the number of subblocks per AWH stored
635 * in the first subblock (we cannot get this directly from the tpr),
636 * so we have to read an AWH block before making the legends.
638 if (awhReader == nullptr)
640 AwhGraphSelection awhGraphSelection = (moreGraphs ? AwhGraphSelection::All : AwhGraphSelection::Pmf);
641 EnergyUnit energyUnit = (kTUnit ? EnergyUnit::KT : EnergyUnit::KJPerMol);
643 gmx::compat::make_unique<AwhReader>(ir.awhParams,
646 energyUnit, BOLTZ*ir.opts.ref_t[0],
650 awhReader->processAwhFrame(*block, frame->t, oenv);
653 while (haveFrame && (timeCheck <= 0));
655 fprintf(stderr, "\n");