2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2019,2020, 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.
37 * \brief Helper code for TPR file access.
39 * \author M. Eric Irrgang <ericirrgang@gmail.com>
40 * \ingroup gmxapi_compat
43 #include "gmxapi/compat/tpr.h"
52 #include "gromacs/mdtypes/inputrec.h"
53 #include "gromacs/topology/topology.h"
54 #include "gromacs/mdtypes/state.h"
55 #include "gromacs/fileio/oenv.h"
56 #include "gromacs/fileio/tpxio.h"
57 #include "gromacs/options/timeunitmanager.h"
58 #include "gromacs/utility/cstringutil.h"
59 #include "gromacs/utility/programcontext.h"
61 #include "gmxapi/gmxapi.h"
62 #include "gmxapi/gmxapicompat.h"
63 #include "gmxapi/compat/mdparams.h"
65 using gmxapi::GmxapiType;
67 namespace gmxapicompat
73 explicit TprContents(const std::string& infile) :
74 irInstance_{ std::make_unique<t_inputrec>() },
75 mtop_{ std::make_unique<gmx_mtop_t>() },
76 state_{ std::make_unique<t_state>() }
78 read_tpx_state(infile.c_str(), irInstance_.get(), state_.get(), mtop_.get());
80 ~TprContents() = default;
81 TprContents(TprContents&& source) noexcept = default;
82 TprContents& operator=(TprContents&&) noexcept = default;
85 * \brief Get a reference to the input record in the TPR file.
87 * Note that this implementation allows different objects to share ownership
88 * of the TprFile and does not provide access restrictions to prevent multiple
89 * code blocks writing to the input record. This should be resolved with a
90 * combination of managed access controlled handles and through better
91 * management of the data structures in the TPR file. I.e. the t_inputrec is
92 * not copyable, moveable, nor default constructable (at least, to produce a
93 * valid record), and it does not necessarily make sense to map the library
94 * data structure to the file data structure (except that we don't have another
95 * way of constructing a complete and valid input record).
97 * \todo We can't play fast and loose with the irInstance for long...
101 t_inputrec& inputRecord() const
107 gmx_mtop_t& molecularTopology() const
113 t_state& state() const
120 // These types are not moveable in GROMACS 2019, so we use unique_ptr as a
121 // moveable wrapper to let TprContents be moveable.
122 std::unique_ptr<t_inputrec> irInstance_;
123 std::unique_ptr<gmx_mtop_t> mtop_;
124 std::unique_ptr<t_state> state_;
127 // Note: This mapping is incomplete. Hopefully we can replace it before more mapping is necessary.
128 // TODO: (#2993) Replace with GROMACS library header resources when available.
129 std::map<std::string, GmxapiType> simulationParameterTypeMap()
132 { "integrator", GmxapiType::STRING },
133 { "tinit", GmxapiType::FLOAT64 },
134 { "dt", GmxapiType::FLOAT64 },
135 { "nsteps", GmxapiType::INT64 },
136 { "init-step", GmxapiType::INT64 },
137 { "simulation-part", GmxapiType::INT64 },
138 { "comm-mode", GmxapiType::STRING },
139 { "nstcomm", GmxapiType::INT64 },
140 { "comm-grps", GmxapiType::NDARRAY }, // Note: we do not have processing for this yet.
141 { "bd-fric", GmxapiType::FLOAT64 },
142 { "ld-seed", GmxapiType::INT64 },
143 { "emtol", GmxapiType::FLOAT64 },
144 { "emstep", GmxapiType::FLOAT64 },
145 { "niter", GmxapiType::INT64 },
146 { "fcstep", GmxapiType::FLOAT64 },
147 { "nstcgsteep", GmxapiType::INT64 },
148 { "nbfgscorr", GmxapiType::INT64 },
149 { "rtpi", GmxapiType::FLOAT64 },
150 { "nstxout", GmxapiType::INT64 },
151 { "nstvout", GmxapiType::INT64 },
152 { "nstfout", GmxapiType::INT64 },
153 { "nstlog", GmxapiType::INT64 },
154 { "nstcalcenergy", GmxapiType::INT64 },
155 { "nstenergy", GmxapiType::INT64 },
156 { "nstxout-compressed", GmxapiType::INT64 },
157 { "compressed-x-precision", GmxapiType::FLOAT64 },
158 { "cutoff-scheme", GmxapiType::STRING },
159 { "nstlist", GmxapiType::INT64 },
160 { "ns-type", GmxapiType::STRING },
161 { "pbc", GmxapiType::STRING },
162 { "periodic-molecules", GmxapiType::BOOL },
169 * TODO: Visitor for predetermined known types.
171 * Development sequence:
174 * 3. template the Visitor setter for compile-time extensibility of type and to prune incompatible types.
175 * 4. switch to Variant type for handling (setter templated on caller input)
176 * 5. switch to Variant type for input as well? (Variant in public API?)
179 std::map<std::string, bool t_inputrec::*> boolParams()
182 { "periodic-molecules", &t_inputrec::bPeriodicMols },
187 std::map<std::string, int t_inputrec::*> int32Params()
190 { "simulation-part", &t_inputrec::simulation_part },
191 { "nstcomm", &t_inputrec::nstcomm },
192 { "niter", &t_inputrec::niter },
193 { "nstcgsteep", &t_inputrec::nstcgsteep },
194 { "nbfgscorr", &t_inputrec::nbfgscorr },
195 { "nstxout", &t_inputrec::nstxout },
196 { "nstvout", &t_inputrec::nstvout },
197 { "nstfout", &t_inputrec::nstfout },
198 { "nstlog", &t_inputrec::nstlog },
199 { "nstcalcenergy", &t_inputrec::nstcalcenergy },
200 { "nstenergy", &t_inputrec::nstenergy },
201 { "nstxout-compressed", &t_inputrec::nstxout_compressed },
202 { "nstlist", &t_inputrec::nstlist },
210 // Provide a helper to disambiguate `real` typed inputrec values.
212 // Primary template returns empty map.
213 template<typename RealT>
214 std::map<std::string, RealT t_inputrec::*> compatibleRealParams()
219 // Specialize for the compile-time typedef of `real` to get the inputrec parameters
220 // compatible with the template parameter whose type is not known until compile time.
222 std::map<std::string, real t_inputrec::*> compatibleRealParams()
225 { "bd-fric", &t_inputrec::bd_fric },
226 { "emtol", &t_inputrec::em_tol },
227 { "emstep", &t_inputrec::em_stepsize },
228 { "fcstep", &t_inputrec::fc_stepsize },
229 { "rtpi", &t_inputrec::rtpi },
230 { "compressed-x-precision", &t_inputrec::x_compression_precision },
238 std::map<std::string, float t_inputrec::*> float32Params()
240 return compatibleRealParams<float>();
243 std::map<std::string, double t_inputrec::*> float64Params()
245 static const std::map<std::string, double t_inputrec::*> explicitDoubles = {
246 { "dt", &t_inputrec::delta_t }, { "tinit", &t_inputrec::init_t },
251 // Initialize map to be returned with any compile-time-only doubles.
252 auto fullMap = compatibleRealParams<double>();
254 // Get the explicitly `double` parameters.
255 for (const auto& item : explicitDoubles)
257 fullMap.emplace(item);
263 std::map<std::string, int64_t t_inputrec::*> int64Params()
266 { "nsteps", &t_inputrec::nsteps },
267 { "init-step", &t_inputrec::init_step },
268 { "ld-seed", &t_inputrec::ld_seed },
275 * \brief Static mapping of parameter names to gmxapi types for GROMACS 2019.
277 * \param name MDP entry name.
278 * \return enumeration value for known parameters.
280 * \throws gmxapi_compat::ValueError for parameters with no mapping.
282 GmxapiType mdParamToType(const std::string& name)
284 const auto staticMap = simulationParameterTypeMap();
285 auto entry = staticMap.find(name);
286 if (entry == staticMap.end())
288 throw ValueError("Named parameter has unknown type mapping.");
290 return entry->second;
295 * \brief Handle / manager for GROMACS molecular computation input parameters.
297 * Interface should be consistent with MDP file entries, but data maps to TPR
298 * file interface. For type safety and simplicity, we don't have generic operator
299 * accessors. Instead, we have templated accessors that throw exceptions when
302 * When MDP input is entirely stored in a key-value tree, this class can be a
303 * simple adapter or wrapper. Until then, we need a manually maintained mapping
304 * of MDP entries to TPR data.
306 * Alternatively, we could update the infrastructure used by list_tpx to provide
307 * more generic output, but our efforts may be better spent in updating the
308 * infrastructure for the key-value tree input system.
310 class GmxMdParamsImpl final
314 * \brief Create an initialized but empty parameters structure.
316 * Parameter keys are set at construction, but all values are empty. This
317 * allows the caller to check for valid parameter names or their types,
318 * while allowing the consuming code to know which parameters were explicitly
321 * To load values from a TPR file, see getMdParams().
325 explicit GmxMdParamsImpl(std::shared_ptr<TprContents> tprContents);
328 * \brief Get the current list of keys.
332 std::vector<std::string> keys() const
334 std::vector<std::string> keyList;
335 for (auto&& entry : int64Params_)
337 keyList.emplace_back(entry.first);
339 for (auto&& entry : intParams_)
341 keyList.emplace_back(entry.first);
343 for (auto&& entry : floatParams_)
345 keyList.emplace_back(entry.first);
347 for (auto&& entry : float64Params_)
349 keyList.emplace_back(entry.first);
355 T extract(const std::string& /* key */) const
358 // should be an APIError
359 throw TypeError("unhandled type");
362 void set(const std::string& key, const int64_t& value)
364 if (int64Params_.find(key) != int64Params_.end())
366 int64Params_[key] = std::make_pair(value, true);
370 auto memberPointer = int64Params().at(key);
371 source_->inputRecord().*memberPointer = value;
374 else if (intParams_.find(key) != intParams_.end())
376 // TODO: check whether value is too large?
377 intParams_[key] = std::make_pair(static_cast<int>(value), true);
381 auto memberPointer = int32Params().at(key);
382 source_->inputRecord().*memberPointer = value;
387 throw KeyError("Named parameter is incompatible with integer type value.");
391 void set(const std::string& key, const double& value)
393 if (float64Params_.find(key) != float64Params_.end())
395 float64Params_[key] = std::make_pair(value, true);
399 auto memberPointer = float64Params().at(key);
400 source_->inputRecord().*memberPointer = value;
403 else if (floatParams_.find(key) != floatParams_.end())
405 // TODO: check whether value is too large?
406 floatParams_[key] = std::make_pair(static_cast<float>(value), true);
410 auto memberPointer = float32Params().at(key);
411 source_->inputRecord().*memberPointer = static_cast<float>(value);
416 throw KeyError("Named parameter is incompatible with floating point type value.");
420 TprReadHandle getSource() const
422 // Note: might return a null handle. Need to decide what that means and how to address it.
423 return TprReadHandle(source_);
427 // Hold the settable parameters and whether or not they have been set.
428 // TODO: update to gmxapi named types?
429 // TODO: update to gmx::compat::optional now that this file is in the GROMACS source.
430 std::map<std::string, std::pair<int64_t, bool>> int64Params_;
431 std::map<std::string, std::pair<int, bool>> intParams_;
432 std::map<std::string, std::pair<float, bool>> floatParams_;
433 std::map<std::string, std::pair<double, bool>> float64Params_;
435 /*! \brief Shared ownership of a pack of TPR data.
437 * This is a non-normative way to retain access to gmxapi resources.
438 * \todo Subscribe to a Context-managed resource.
440 std::shared_ptr<TprContents> source_;
443 void setParam(gmxapicompat::GmxMdParams* params, const std::string& name, double value)
445 assert(params != nullptr);
446 assert(params->params_ != nullptr);
447 params->params_->set(name, value);
450 void setParam(gmxapicompat::GmxMdParams* params, const std::string& name, int64_t value)
452 assert(params != nullptr);
453 assert(params->params_ != nullptr);
454 params->params_->set(name, value);
457 template<typename ParamsContainerT, typename Mapping>
458 static void updateParamsContainer(ParamsContainerT* params, const TprContents& source, const Mapping& map)
460 for (const auto& definition : map)
462 const auto& key = definition.first;
463 auto memberPointer = definition.second;
464 auto& irInstance = source.inputRecord();
465 auto fileValue = irInstance.*memberPointer;
466 (*params)[key] = std::make_pair(fileValue, true);
471 * \brief A GmxMdParams implementation that depends on TPR files.
475 GmxMdParamsImpl::GmxMdParamsImpl(std::shared_ptr<gmxapicompat::TprContents> tprContents) :
476 source_{ std::move(tprContents) }
480 updateParamsContainer(&int64Params_, *source_, int64Params());
481 updateParamsContainer(&intParams_, *source_, int32Params());
482 updateParamsContainer(&floatParams_, *source_, float32Params());
483 updateParamsContainer(&float64Params_, *source_, float64Params());
487 GmxMdParamsImpl::GmxMdParamsImpl() : GmxMdParamsImpl(nullptr) {}
490 int GmxMdParamsImpl::extract<int>(const std::string& key) const
492 const auto& params = intParams_;
493 const auto& entry = params.find(key);
494 if (entry == params.cend())
496 throw KeyError("Parameter of the requested name and type not defined.");
498 else if (!entry->second.second)
500 // TODO: handle invalid and unset parameters differently.
501 throw KeyError("Parameter of the requested name not set.");
505 return entry->second.first;
510 int64_t GmxMdParamsImpl::extract<int64_t>(const std::string& key) const
512 const auto& params = int64Params_;
513 const auto& entry = params.find(key);
514 if (entry == params.cend())
516 throw KeyError("Parameter of the requested name and type not defined.");
518 else if (!entry->second.second)
520 // TODO: handle invalid and unset parameters differently.
521 throw KeyError("Parameter of the requested name not set.");
525 return entry->second.first;
529 float GmxMdParamsImpl::extract<float>(const std::string& key) const
531 const auto& params = floatParams_;
532 const auto& entry = params.find(key);
533 if (entry == params.cend())
535 throw KeyError("Parameter of the requested name and type not defined.");
537 else if (!entry->second.second)
539 // TODO: handle invalid and unset parameters differently.
540 throw KeyError("Parameter of the requested name not set.");
544 return entry->second.first;
548 double GmxMdParamsImpl::extract<double>(const std::string& key) const
550 const auto& params = float64Params_;
551 const auto& entry = params.find(key);
552 if (entry == params.cend())
554 throw KeyError("Parameter of the requested name and type not defined.");
556 else if (!entry->second.second)
558 // TODO: handle invalid and unset parameters differently.
559 throw KeyError("Parameter of the requested name not set.");
563 return entry->second.first;
568 int extractParam(const GmxMdParams& params, const std::string& name, int /*unused*/)
570 assert(params.params_);
571 return params.params_->extract<int>(name);
574 int64_t extractParam(const GmxMdParams& params, const std::string& name, int64_t /*unused*/)
576 assert(params.params_);
578 // Allow fetching both known integer types.
581 value = params.params_->extract<int>(name);
583 catch (const KeyError& error)
585 // If not found as a regular int, check for int64.
588 value = params.params_->extract<int64_t>(name);
590 catch (const KeyError& error64)
592 throw KeyError("Parameter of the requested name not set.");
595 // Any other exceptions propagate out.
599 float extractParam(const GmxMdParams& params, const std::string& name, float /*unused*/)
601 assert(params.params_);
602 return params.params_->extract<float>(name);
605 double extractParam(const GmxMdParams& params, const std::string& name, double /*unused*/)
607 assert(params.params_);
609 // Allow fetching both single and double precision.
612 value = params.params_->extract<double>(name);
614 catch (const KeyError& errorDouble)
616 // If not found as a double precision value, check for single-precision.
619 value = params.params_->extract<float>(name);
621 catch (const KeyError& errorFloat)
623 throw KeyError("Parameter of the requested name not set.");
626 // Any other exceptions propagate out.
630 std::vector<std::string> keys(const GmxMdParams& params)
632 return params.params_->keys();
636 std::unique_ptr<TprReadHandle> readTprFile(const std::string& filename)
638 auto tprfile = gmxapicompat::TprContents(filename);
639 auto handle = std::make_unique<gmxapicompat::TprReadHandle>(std::move(tprfile));
643 std::unique_ptr<GmxMdParams> getMdParams(const TprReadHandle& handle)
645 auto tprfile = handle.get();
646 // TODO: convert to exception / decide whether null handles are allowed.
648 auto params_impl = std::make_unique<GmxMdParamsImpl>(tprfile);
649 auto params = std::make_unique<GmxMdParams>(std::move(params_impl));
653 std::unique_ptr<TopologySource> getTopologySource(const TprReadHandle& handle)
655 auto source = std::make_unique<TopologySource>();
656 source->tprFile_ = handle.get();
660 std::unique_ptr<SimulationState> getSimulationState(const TprReadHandle& handle)
662 auto source = std::make_unique<SimulationState>();
663 source->tprFile_ = handle.get();
667 std::unique_ptr<StructureSource> getStructureSource(const TprReadHandle& handle)
669 auto source = std::make_unique<StructureSource>();
670 source->tprFile_ = handle.get();
674 TprReadHandle::TprReadHandle(std::shared_ptr<TprContents> tprFile) :
675 tprContents_{ std::move(tprFile) }
679 TprReadHandle getSourceFileHandle(const GmxMdParams& params)
681 return params.params_->getSource();
684 void writeTprFile(const std::string& filename,
685 const GmxMdParams& params,
686 const StructureSource& structure,
687 const SimulationState& state,
688 const TopologySource& topology)
690 assert(params.params_);
691 // The only way we can check for consistent input right now is to make sure
692 // it all comes from the same file.
693 if (structure.tprFile_.get() != state.tprFile_.get()
694 || state.tprFile_.get() != topology.tprFile_.get()
695 || topology.tprFile_.get() != params.params_->getSource().get().get()
696 || params.params_->getSource().get().get() != structure.tprFile_.get())
699 "writeTprFile does not yet know how to reconcile data from different TPR file "
703 const auto tprFileHandle = params.params_->getSource();
704 const auto tprFile = tprFileHandle.get();
706 const auto& inputRecord = tprFile->inputRecord();
707 const auto& writeState = tprFile->state();
708 const auto& writeTopology = tprFile->molecularTopology();
709 write_tpx_state(filename.c_str(), &inputRecord, &writeState, &writeTopology);
712 TprReadHandle::TprReadHandle(TprContents&& tprFile) :
713 TprReadHandle{ std::make_shared<TprContents>(std::move(tprFile)) }
717 std::shared_ptr<TprContents> TprReadHandle::get() const
722 // defaulted here to delay definition until after member types are defined.
723 TprReadHandle::~TprReadHandle() = default;
725 GmxMdParams::~GmxMdParams() = default;
727 GmxMdParams::GmxMdParams() : params_{ std::make_unique<GmxMdParamsImpl>() } {}
729 GmxMdParams::GmxMdParams(GmxMdParams&&) noexcept = default;
731 GmxMdParams& GmxMdParams::operator=(GmxMdParams&&) noexcept = default;
733 GmxMdParams::GmxMdParams(std::unique_ptr<GmxMdParamsImpl>&& impl)
735 // We use swap instead of move construction so that we don't have
736 // to worry about the restrictions on Deleters.
737 // Ref: https://en.cppreference.com/w/cpp/memory/unique_ptr/unique_ptr
741 // maybe this should return a handle to the new file?
742 bool copy_tprfile(const gmxapicompat::TprReadHandle& input, const std::string& outFile)
748 gmxapicompat::writeTprFile(
749 outFile, *gmxapicompat::getMdParams(input), *gmxapicompat::getStructureSource(input),
750 *gmxapicompat::getSimulationState(input), *gmxapicompat::getTopologySource(input));
754 bool rewrite_tprfile(const std::string& inFile, const std::string& outFile, double endTime)
756 auto success = false;
758 const char* top_fn = inFile.c_str();
760 t_inputrec irInstance;
763 read_tpx_state(top_fn, &irInstance, &state, &mtop);
765 /* set program name, command line, and default values for output options */
766 gmx_output_env_t* oenv;
767 gmx::TimeUnit timeUnit = gmx::TimeUnit::Default;
768 bool bView{ false }; // argument that says we don't want to view graphs.
769 output_env_init(&oenv, gmx::getProgramContext(), timeUnit, bView, XvgFormat::Xmgrace, 0);
771 double run_t = irInstance.init_step * irInstance.delta_t + irInstance.init_t;
773 irInstance.nsteps = lround((endTime - run_t) / irInstance.delta_t);
775 write_tpx_state(outFile.c_str(), &irInstance, &state, &mtop);
781 } // end namespace gmxapicompat