api/gmxapi/cpp/tpr.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2019,2020,2021, 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 /*! \file
  37  * \brief Helper code for TPR file access.
  38  *
  39  * \author M. Eric Irrgang <ericirrgang@gmail.com>
  40  * \ingroup gmxapi_compat
  41  */
  42
  43 #include "gmxapi/compat/tpr.h"
  44
  45 #include <cassert>
  46
  47 #include <map>
  48 #include <memory>
  49 #include <string>
  50 #include <vector>
  51
  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"
  60
  61 #include "gmxapi/gmxapi.h"
  62 #include "gmxapi/gmxapicompat.h"
  63 #include "gmxapi/compat/mdparams.h"
  64
  65 using gmxapi::GmxapiType;
  66
  67 namespace gmxapicompat
  68 {
  69
  70 class TprContents
  71 {
  72 public:
  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>() }
  77     {
  78         read_tpx_state(infile.c_str(), irInstance_.get(), state_.get(), mtop_.get());
  79     }
  80     ~TprContents()                             = default;
  81     TprContents(TprContents&& source) noexcept = default;
  82     TprContents& operator=(TprContents&&) noexcept = default;
  83
  84     /*!
  85      * \brief Get a reference to the input record in the TPR file.
  86      *
  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).
  96      *
  97      * \todo We can't play fast and loose with the irInstance for long...
  98      *
  99      * \return
 100      */
 101     t_inputrec& inputRecord() const
 102     {
 103         assert(irInstance_);
 104         return *irInstance_;
 105     }
 106
 107     gmx_mtop_t& molecularTopology() const
 108     {
 109         assert(mtop_);
 110         return *mtop_;
 111     }
 112
 113     t_state& state() const
 114     {
 115         assert(state_);
 116         return *state_;
 117     }
 118
 119 private:
 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_;
 125 };
 126
 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()
 130 {
 131     return {
 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 },
 163         //            TBD
 164
 165     };
 166 }
 167
 168 /*
 169  * TODO: Visitor for predetermined known types.
 170  *
 171  * Development sequence:
 172  * 1. map pointers
 173  * 2. map setters ()
 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?)
 177  */
 178
 179 std::map<std::string, bool t_inputrec::*> boolParams()
 180 {
 181     return {
 182         { "periodic-molecules", &t_inputrec::bPeriodicMols },
 183         //            ...
 184     };
 185 }
 186
 187 std::map<std::string, int t_inputrec::*> int32Params()
 188 {
 189     return {
 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 },
 203         //            ...
 204     };
 205 }
 206
 207 namespace
 208 {
 209
 210 // Provide a helper to disambiguate `real` typed inputrec values.
 211
 212 // Primary template returns empty map.
 213 template<typename RealT>
 214 std::map<std::string, RealT t_inputrec::*> compatibleRealParams()
 215 {
 216     return {};
 217 }
 218
 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.
 221 template<>
 222 std::map<std::string, real t_inputrec::*> compatibleRealParams()
 223 {
 224     return {
 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 },
 231         //            ...
 232
 233     };
 234 }
 235
 236 } // namespace
 237
 238 std::map<std::string, float t_inputrec::*> float32Params()
 239 {
 240     return compatibleRealParams<float>();
 241 }
 242
 243 std::map<std::string, double t_inputrec::*> float64Params()
 244 {
 245     static const std::map<std::string, double t_inputrec::*> explicitDoubles = {
 246         { "dt", &t_inputrec::delta_t }, { "tinit", &t_inputrec::init_t },
 247         //            ...
 248
 249     };
 250
 251     // Initialize map to be returned with any compile-time-only doubles.
 252     auto fullMap = compatibleRealParams<double>();
 253
 254     // Get the explicitly `double` parameters.
 255     for (const auto& item : explicitDoubles)
 256     {
 257         fullMap.emplace(item);
 258     }
 259
 260     return fullMap;
 261 }
 262
 263 std::map<std::string, int64_t t_inputrec::*> int64Params()
 264 {
 265     return {
 266         { "nsteps", &t_inputrec::nsteps },
 267         { "init-step", &t_inputrec::init_step },
 268         { "ld-seed", &t_inputrec::ld_seed },
 269         //            ...
 270
 271     };
 272 }
 273
 274 /*!
 275  * \brief Static mapping of parameter names to gmxapi types for GROMACS 2019.
 276  *
 277  * \param name MDP entry name.
 278  * \return enumeration value for known parameters.
 279  *
 280  * \throws gmxapi_compat::ValueError for parameters with no mapping.
 281  */
 282 GmxapiType mdParamToType(const std::string& name)
 283 {
 284     const auto staticMap = simulationParameterTypeMap();
 285     auto       entry     = staticMap.find(name);
 286     if (entry == staticMap.end())
 287     {
 288         throw ValueError("Named parameter has unknown type mapping.");
 289     }
 290     return entry->second;
 291 };
 292
 293
 294 /*!
 295  * \brief Handle / manager for GROMACS molecular computation input parameters.
 296  *
 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
 300  * there is trouble.
 301  *
 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.
 305  *
 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.
 309  */
 310 class GmxMdParamsImpl final
 311 {
 312 public:
 313     /*!
 314      * \brief Create an initialized but empty parameters structure.
 315      *
 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
 319      * set by the caller.
 320      *
 321      * To load values from a TPR file, see getMdParams().
 322      */
 323     GmxMdParamsImpl();
 324
 325     explicit GmxMdParamsImpl(std::shared_ptr<TprContents> tprContents);
 326
 327     /*!
 328      * \brief Get the current list of keys.
 329      *
 330      * \return
 331      */
 332     std::vector<std::string> keys() const
 333     {
 334         std::vector<std::string> keyList;
 335         for (auto&& entry : int64Params_)
 336         {
 337             keyList.emplace_back(entry.first);
 338         }
 339         for (auto&& entry : intParams_)
 340         {
 341             keyList.emplace_back(entry.first);
 342         }
 343         for (auto&& entry : floatParams_)
 344         {
 345             keyList.emplace_back(entry.first);
 346         }
 347         for (auto&& entry : float64Params_)
 348         {
 349             keyList.emplace_back(entry.first);
 350         }
 351         return keyList;
 352     };
 353
 354     template<typename T>
 355     T extract(const std::string& /* key */) const
 356     {
 357         auto value = T();
 358         // should be an APIError
 359         throw TypeError("unhandled type");
 360     }
 361
 362     void set(const std::string& key, const int64_t& value)
 363     {
 364         if (int64Params_.find(key) != int64Params_.end())
 365         {
 366             int64Params_[key] = std::make_pair(value, true);
 367
 368             if (source_)
 369             {
 370                 auto memberPointer                    = int64Params().at(key);
 371                 source_->inputRecord().*memberPointer = value;
 372             }
 373         }
 374         else if (intParams_.find(key) != intParams_.end())
 375         {
 376             // TODO: check whether value is too large?
 377             intParams_[key] = std::make_pair(static_cast<int>(value), true);
 378
 379             if (source_)
 380             {
 381                 auto memberPointer                    = int32Params().at(key);
 382                 source_->inputRecord().*memberPointer = value;
 383             }
 384         }
 385         else
 386         {
 387             throw KeyError("Named parameter is incompatible with integer type value.");
 388         }
 389     };
 390
 391     void set(const std::string& key, const double& value)
 392     {
 393         if (float64Params_.find(key) != float64Params_.end())
 394         {
 395             float64Params_[key] = std::make_pair(value, true);
 396
 397             if (source_)
 398             {
 399                 auto memberPointer                    = float64Params().at(key);
 400                 source_->inputRecord().*memberPointer = value;
 401             }
 402         }
 403         else if (floatParams_.find(key) != floatParams_.end())
 404         {
 405             // TODO: check whether value is too large?
 406             floatParams_[key] = std::make_pair(static_cast<float>(value), true);
 407
 408             if (source_)
 409             {
 410                 auto memberPointer                    = float32Params().at(key);
 411                 source_->inputRecord().*memberPointer = static_cast<float>(value);
 412             }
 413         }
 414         else
 415         {
 416             throw KeyError("Named parameter is incompatible with floating point type value.");
 417         }
 418     };
 419
 420     TprReadHandle getSource() const
 421     {
 422         // Note: might return a null handle. Need to decide what that means and how to address it.
 423         return TprReadHandle(source_);
 424     }
 425
 426 private:
 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_;
 434
 435     /*! \brief Shared ownership of a pack of TPR data.
 436      *
 437      * This is a non-normative way to retain access to gmxapi resources.
 438      * \todo Subscribe to a Context-managed resource.
 439      */
 440     std::shared_ptr<TprContents> source_;
 441 };
 442
 443 void setParam(gmxapicompat::GmxMdParams* params, const std::string& name, double value)
 444 {
 445     assert(params != nullptr);
 446     assert(params->params_ != nullptr);
 447     params->params_->set(name, value);
 448 }
 449
 450 void setParam(gmxapicompat::GmxMdParams* params, const std::string& name, int64_t value)
 451 {
 452     assert(params != nullptr);
 453     assert(params->params_ != nullptr);
 454     params->params_->set(name, value);
 455 }
 456
 457 template<typename ParamsContainerT, typename Mapping>
 458 static void updateParamsContainer(ParamsContainerT* params, const TprContents& source, const Mapping& map)
 459 {
 460     for (const auto& definition : map)
 461     {
 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);
 467     }
 468 }
 469
 470 /*!
 471  * \brief A GmxMdParams implementation that depends on TPR files.
 472  *
 473  * \param tprContents
 474  */
 475 GmxMdParamsImpl::GmxMdParamsImpl(std::shared_ptr<gmxapicompat::TprContents> tprContents) :
 476     source_{ std::move(tprContents) }
 477 {
 478     if (source_)
 479     {
 480         updateParamsContainer(&int64Params_, *source_, int64Params());
 481         updateParamsContainer(&intParams_, *source_, int32Params());
 482         updateParamsContainer(&floatParams_, *source_, float32Params());
 483         updateParamsContainer(&float64Params_, *source_, float64Params());
 484     }
 485 }
 486
 487 GmxMdParamsImpl::GmxMdParamsImpl() : GmxMdParamsImpl(nullptr) {}
 488
 489 template<>
 490 int GmxMdParamsImpl::extract<int>(const std::string& key) const
 491 {
 492     const auto& params = intParams_;
 493     const auto& entry  = params.find(key);
 494     if (entry == params.cend())
 495     {
 496         throw KeyError("Parameter of the requested name and type not defined.");
 497     }
 498     else if (!entry->second.second)
 499     {
 500         // TODO: handle invalid and unset parameters differently.
 501         throw KeyError("Parameter of the requested name not set.");
 502     }
 503     else
 504     {
 505         return entry->second.first;
 506     }
 507 }
 508
 509 template<>
 510 int64_t GmxMdParamsImpl::extract<int64_t>(const std::string& key) const
 511 {
 512     const auto& params = int64Params_;
 513     const auto& entry  = params.find(key);
 514     if (entry == params.cend())
 515     {
 516         throw KeyError("Parameter of the requested name and type not defined.");
 517     }
 518     else if (!entry->second.second)
 519     {
 520         // TODO: handle invalid and unset parameters differently.
 521         throw KeyError("Parameter of the requested name not set.");
 522     }
 523     else
 524     {
 525         return entry->second.first;
 526     }
 527 }
 528 template<>
 529 float GmxMdParamsImpl::extract<float>(const std::string& key) const
 530 {
 531     const auto& params = floatParams_;
 532     const auto& entry  = params.find(key);
 533     if (entry == params.cend())
 534     {
 535         throw KeyError("Parameter of the requested name and type not defined.");
 536     }
 537     else if (!entry->second.second)
 538     {
 539         // TODO: handle invalid and unset parameters differently.
 540         throw KeyError("Parameter of the requested name not set.");
 541     }
 542     else
 543     {
 544         return entry->second.first;
 545     }
 546 }
 547 template<>
 548 double GmxMdParamsImpl::extract<double>(const std::string& key) const
 549 {
 550     const auto& params = float64Params_;
 551     const auto& entry  = params.find(key);
 552     if (entry == params.cend())
 553     {
 554         throw KeyError("Parameter of the requested name and type not defined.");
 555     }
 556     else if (!entry->second.second)
 557     {
 558         // TODO: handle invalid and unset parameters differently.
 559         throw KeyError("Parameter of the requested name not set.");
 560     }
 561     else
 562     {
 563         return entry->second.first;
 564     }
 565 }
 566
 567
 568 int extractParam(const GmxMdParams& params, const std::string& name, int /*unused*/)
 569 {
 570     assert(params.params_);
 571     return params.params_->extract<int>(name);
 572 }
 573
 574 int64_t extractParam(const GmxMdParams& params, const std::string& name, int64_t /*unused*/)
 575 {
 576     assert(params.params_);
 577     int64_t value{};
 578     // Allow fetching both known integer types.
 579     try
 580     {
 581         value = params.params_->extract<int>(name);
 582     }
 583     catch (const KeyError& error)
 584     {
 585         // If not found as a regular int, check for int64.
 586         try
 587         {
 588             value = params.params_->extract<int64_t>(name);
 589         }
 590         catch (const KeyError& error64)
 591         {
 592             throw KeyError("Parameter of the requested name not set.");
 593         }
 594     }
 595     // Any other exceptions propagate out.
 596     return value;
 597 }
 598
 599 float extractParam(const GmxMdParams& params, const std::string& name, float /*unused*/)
 600 {
 601     assert(params.params_);
 602     return params.params_->extract<float>(name);
 603 }
 604
 605 double extractParam(const GmxMdParams& params, const std::string& name, double /*unused*/)
 606 {
 607     assert(params.params_);
 608     double value{};
 609     // Allow fetching both single and double precision.
 610     try
 611     {
 612         value = params.params_->extract<double>(name);
 613     }
 614     catch (const KeyError& errorDouble)
 615     {
 616         // If not found as a double precision value, check for single-precision.
 617         try
 618         {
 619             value = params.params_->extract<float>(name);
 620         }
 621         catch (const KeyError& errorFloat)
 622         {
 623             throw KeyError("Parameter of the requested name not set.");
 624         }
 625     }
 626     // Any other exceptions propagate out.
 627     return value;
 628 }
 629
 630 std::vector<std::string> keys(const GmxMdParams& params)
 631 {
 632     return params.params_->keys();
 633 }
 634
 635
 636 std::unique_ptr<TprReadHandle> readTprFile(const std::string& filename)
 637 {
 638     auto tprfile = gmxapicompat::TprContents(filename);
 639     auto handle  = std::make_unique<gmxapicompat::TprReadHandle>(std::move(tprfile));
 640     return handle;
 641 }
 642
 643 std::unique_ptr<GmxMdParams> getMdParams(const TprReadHandle& handle)
 644 {
 645     auto tprfile = handle.get();
 646     // TODO: convert to exception / decide whether null handles are allowed.
 647     assert(tprfile);
 648     auto params_impl = std::make_unique<GmxMdParamsImpl>(tprfile);
 649     auto params      = std::make_unique<GmxMdParams>(std::move(params_impl));
 650     return params;
 651 }
 652
 653 std::unique_ptr<TopologySource> getTopologySource(const TprReadHandle& handle)
 654 {
 655     auto source      = std::make_unique<TopologySource>();
 656     source->tprFile_ = handle.get();
 657     return source;
 658 }
 659
 660 std::unique_ptr<SimulationState> getSimulationState(const TprReadHandle& handle)
 661 {
 662     auto source      = std::make_unique<SimulationState>();
 663     source->tprFile_ = handle.get();
 664     return source;
 665 }
 666
 667 std::unique_ptr<StructureSource> getStructureSource(const TprReadHandle& handle)
 668 {
 669     auto source      = std::make_unique<StructureSource>();
 670     source->tprFile_ = handle.get();
 671     return source;
 672 }
 673
 674 TprReadHandle::TprReadHandle(std::shared_ptr<TprContents> tprFile) :
 675     tprContents_{ std::move(tprFile) }
 676 {
 677 }
 678
 679 TprReadHandle getSourceFileHandle(const GmxMdParams& params)
 680 {
 681     return params.params_->getSource();
 682 }
 683
 684 void writeTprFile(const std::string&     filename,
 685                   const GmxMdParams&     params,
 686                   const StructureSource& structure,
 687                   const SimulationState& state,
 688                   const TopologySource&  topology)
 689 {
 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())
 697     {
 698         throw ValueError(
 699                 "writeTprFile does not yet know how to reconcile data from different TPR file "
 700                 "sources.");
 701     }
 702
 703     const auto tprFileHandle = params.params_->getSource();
 704     const auto tprFile       = tprFileHandle.get();
 705     assert(tprFile);
 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);
 710 }
 711
 712 TprReadHandle::TprReadHandle(TprContents&& tprFile) :
 713     TprReadHandle{ std::make_shared<TprContents>(std::move(tprFile)) }
 714 {
 715 }
 716
 717 std::shared_ptr<TprContents> TprReadHandle::get() const
 718 {
 719     return tprContents_;
 720 }
 721
 722 // defaulted here to delay definition until after member types are defined.
 723 TprReadHandle::~TprReadHandle() = default;
 724
 725 GmxMdParams::~GmxMdParams() = default;
 726
 727 GmxMdParams::GmxMdParams() : params_{ std::make_unique<GmxMdParamsImpl>() } {}
 728
 729 GmxMdParams::GmxMdParams(GmxMdParams&&) noexcept = default;
 730
 731 GmxMdParams& GmxMdParams::operator=(GmxMdParams&&) noexcept = default;
 732
 733 GmxMdParams::GmxMdParams(std::unique_ptr<GmxMdParamsImpl>&& impl)
 734 {
 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
 738     params_.swap(impl);
 739 };
 740
 741 // maybe this should return a handle to the new file?
 742 bool copy_tprfile(const gmxapicompat::TprReadHandle& input, const std::string& outFile)
 743 {
 744     if (!input.get())
 745     {
 746         return false;
 747     }
 748     gmxapicompat::writeTprFile(outFile,
 749                                *gmxapicompat::getMdParams(input),
 750                                *gmxapicompat::getStructureSource(input),
 751                                *gmxapicompat::getSimulationState(input),
 752                                *gmxapicompat::getTopologySource(input));
 753     return true;
 754 }
 755
 756 bool rewrite_tprfile(const std::string& inFile, const std::string& outFile, double endTime)
 757 {
 758     auto success = false;
 759
 760     const char* top_fn = inFile.c_str();
 761
 762     t_inputrec irInstance;
 763     gmx_mtop_t mtop;
 764     t_state    state;
 765     read_tpx_state(top_fn, &irInstance, &state, &mtop);
 766
 767     /* set program name, command line, and default values for output options */
 768     gmx_output_env_t* oenv;
 769     gmx::TimeUnit     timeUnit = gmx::TimeUnit::Default;
 770     bool              bView{ false }; // argument that says we don't want to view graphs.
 771     output_env_init(&oenv, gmx::getProgramContext(), timeUnit, bView, XvgFormat::Xmgrace, 0);
 772
 773     double run_t = irInstance.init_step * irInstance.delta_t + irInstance.init_t;
 774
 775     irInstance.nsteps = lround((endTime - run_t) / irInstance.delta_t);
 776
 777     write_tpx_state(outFile.c_str(), &irInstance, &state, mtop);
 778
 779     success = true;
 780     return success;
 781 }
 782
 783 } // end namespace gmxapicompat