2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2010,2011,2012,2013,2014,2015,2016,2017,2018,2019, 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 * Implements gmx::TrajectoryAnalysisRunnerCommon.
39 * \author Teemu Murtola <teemu.murtola@gmail.com>
40 * \ingroup module_trajectoryanalysis
44 #include "runnercommon.h"
51 #include "gromacs/fileio/oenv.h"
52 #include "gromacs/fileio/timecontrol.h"
53 #include "gromacs/fileio/trxio.h"
54 #include "gromacs/math/vec.h"
55 #include "gromacs/options/basicoptions.h"
56 #include "gromacs/options/filenameoption.h"
57 #include "gromacs/options/ioptionscontainer.h"
58 #include "gromacs/pbcutil/rmpbc.h"
59 #include "gromacs/selection/selection.h"
60 #include "gromacs/selection/selectioncollection.h"
61 #include "gromacs/selection/selectionoption.h"
62 #include "gromacs/selection/selectionoptionbehavior.h"
63 #include "gromacs/topology/topology.h"
64 #include "gromacs/trajectory/trajectoryframe.h"
65 #include "gromacs/trajectoryanalysis/analysissettings.h"
66 #include "gromacs/trajectoryanalysis/topologyinformation.h"
67 #include "gromacs/utility/cstringutil.h"
68 #include "gromacs/utility/exceptions.h"
69 #include "gromacs/utility/gmxassert.h"
70 #include "gromacs/utility/programcontext.h"
71 #include "gromacs/utility/smalloc.h"
72 #include "gromacs/utility/stringutil.h"
74 #include "analysissettings_impl.h"
79 class TrajectoryAnalysisRunnerCommon::Impl : public ITopologyProvider
82 explicit Impl(TrajectoryAnalysisSettings *settings);
85 bool hasTrajectory() const { return !trjfile_.empty(); }
87 void initTopology(bool required);
88 void initFirstFrame();
89 void initFrameIndexGroup();
90 void finishTrajectory();
92 // From ITopologyProvider
93 gmx_mtop_t *getTopology(bool required) override
95 initTopology(required);
96 return topInfo_.mtop_.get();
98 int getAtomCount() override
100 if (!topInfo_.hasTopology())
102 if (trajectoryGroup_.isValid())
104 GMX_THROW(InconsistentInputError("-fgroup is only supported when -s is also specified"));
106 // Read the first frame if we don't know the maximum number of
114 TrajectoryAnalysisSettings &settings_;
115 TopologyInformation topInfo_;
117 //! Name of the trajectory file (empty if not provided).
118 std::string trjfile_;
119 //! Name of the topology file (empty if no topology provided).
120 std::string topfile_;
121 Selection trajectoryGroup_;
130 //! The current frame, or \p NULL if no frame loaded yet.
133 //! Used to store the status variable from read_first_frame().
134 t_trxstatus *status_;
135 gmx_output_env_t *oenv_;
139 TrajectoryAnalysisRunnerCommon::Impl::Impl(TrajectoryAnalysisSettings *settings)
140 : settings_(*settings),
141 startTime_(0.0), endTime_(0.0), deltaTime_(0.0),
142 bStartTimeSet_(false), bEndTimeSet_(false), bDeltaTimeSet_(false),
143 bTrajOpen_(false), fr(nullptr), gpbc_(nullptr), status_(nullptr), oenv_(nullptr)
148 TrajectoryAnalysisRunnerCommon::Impl::~Impl()
153 // There doesn't seem to be a function for freeing frame data
160 if (oenv_ != nullptr)
162 output_env_done(oenv_);
167 TrajectoryAnalysisRunnerCommon::Impl::initTopology(bool required)
169 // Return immediately if the topology has already been loaded.
170 if (topInfo_.hasTopology())
175 if (required && topfile_.empty())
177 GMX_THROW(InconsistentInputError("No topology provided, but one is required for analysis"));
180 // Load the topology if requested.
181 if (!topfile_.empty())
183 topInfo_.fillFromInputFile(topfile_);
185 && !settings_.hasFlag(TrajectoryAnalysisSettings::efUseTopX))
187 topInfo_.xtop_.clear();
190 && !settings_.hasFlag(TrajectoryAnalysisSettings::efUseTopV))
192 topInfo_.vtop_.clear();
198 TrajectoryAnalysisRunnerCommon::Impl::initFirstFrame()
200 // Return if we have already initialized the trajectory.
205 time_unit_t time_unit
206 = static_cast<time_unit_t>(settings_.timeUnit() + 1); // NOLINT(bugprone-misplaced-widening-cast)
207 output_env_init(&oenv_, getProgramContext(), time_unit, FALSE, exvgNONE, 0);
209 int frflags = settings_.frflags();
210 frflags |= TRX_NEED_X;
216 if (!read_first_frame(oenv_, &status_, trjfile_.c_str(), fr, frflags))
218 GMX_THROW(FileIOError("Could not read coordinates from trajectory"));
222 if (topInfo_.hasTopology())
224 const int topologyAtomCount = topInfo_.mtop()->natoms;
225 if (fr->natoms > topologyAtomCount)
227 const std::string message
228 = formatString("Trajectory (%d atoms) does not match topology (%d atoms)",
229 fr->natoms, topologyAtomCount);
230 GMX_THROW(InconsistentInputError(message));
236 // Prepare a frame from topology information.
237 if (frflags & (TRX_NEED_F))
239 GMX_THROW(InvalidInputError("Forces cannot be read from a topology"));
241 fr->natoms = topInfo_.mtop()->natoms;
243 snew(fr->x, fr->natoms);
244 memcpy(fr->x, topInfo_.xtop_.data(),
245 sizeof(*fr->x) * fr->natoms);
246 if (frflags & (TRX_NEED_V))
248 if (topInfo_.vtop_.empty())
250 GMX_THROW(InvalidInputError("Velocities were required, but could not be read from the topology file"));
253 snew(fr->v, fr->natoms);
254 memcpy(fr->v, topInfo_.vtop_.data(),
255 sizeof(*fr->v) * fr->natoms);
258 copy_mat(topInfo_.boxtop_, fr->box);
261 set_trxframe_ePBC(fr, topInfo_.ePBC());
262 if (topInfo_.hasTopology() && settings_.hasRmPBC())
264 gpbc_ = gmx_rmpbc_init(topInfo_);
269 TrajectoryAnalysisRunnerCommon::Impl::initFrameIndexGroup()
271 if (!trajectoryGroup_.isValid())
275 GMX_RELEASE_ASSERT(bTrajOpen_,
276 "Trajectory index only makes sense with a real trajectory");
277 if (trajectoryGroup_.atomCount() != fr->natoms)
279 const std::string message = formatString(
280 "Selection specified with -fgroup has %d atoms, but "
281 "the trajectory (-f) has %d atoms.",
282 trajectoryGroup_.atomCount(), fr->natoms);
283 GMX_THROW(InconsistentInputError(message));
286 snew(fr->index, trajectoryGroup_.atomCount());
287 std::copy(trajectoryGroup_.atomIndices().begin(),
288 trajectoryGroup_.atomIndices().end(),
293 TrajectoryAnalysisRunnerCommon::Impl::finishTrajectory()
300 if (gpbc_ != nullptr)
302 gmx_rmpbc_done(gpbc_);
307 /*********************************************************************
308 * TrajectoryAnalysisRunnerCommon
311 TrajectoryAnalysisRunnerCommon::TrajectoryAnalysisRunnerCommon(
312 TrajectoryAnalysisSettings *settings)
313 : impl_(new Impl(settings))
318 TrajectoryAnalysisRunnerCommon::~TrajectoryAnalysisRunnerCommon()
324 TrajectoryAnalysisRunnerCommon::topologyProvider()
331 TrajectoryAnalysisRunnerCommon::initOptions(IOptionsContainer *options,
332 TimeUnitBehavior *timeUnitBehavior)
334 TrajectoryAnalysisSettings &settings = impl_->settings_;
336 // Add common file name arguments.
337 options->addOption(FileNameOption("f")
338 .filetype(eftTrajectory).inputFile()
339 .store(&impl_->trjfile_)
340 .defaultBasename("traj")
341 .description("Input trajectory or single configuration"));
342 options->addOption(FileNameOption("s")
343 .filetype(eftTopology).inputFile()
344 .store(&impl_->topfile_)
345 .defaultBasename("topol")
346 .description("Input structure"));
348 // Add options for trajectory time control.
349 options->addOption(DoubleOption("b")
350 .store(&impl_->startTime_).storeIsSet(&impl_->bStartTimeSet_)
352 .description("First frame (%t) to read from trajectory"));
353 options->addOption(DoubleOption("e")
354 .store(&impl_->endTime_).storeIsSet(&impl_->bEndTimeSet_)
356 .description("Last frame (%t) to read from trajectory"));
357 options->addOption(DoubleOption("dt")
358 .store(&impl_->deltaTime_).storeIsSet(&impl_->bDeltaTimeSet_)
360 .description("Only use frame if t MOD dt == first time (%t)"));
362 // Add time unit option.
363 timeUnitBehavior->setTimeUnitFromEnvironment();
364 timeUnitBehavior->addTimeUnitOption(options, "tu");
365 timeUnitBehavior->setTimeUnitStore(&impl_->settings_.impl_->timeUnit);
367 options->addOption(SelectionOption("fgroup")
368 .store(&impl_->trajectoryGroup_)
369 .onlySortedAtoms().onlyStatic()
370 .description("Atoms stored in the trajectory file "
371 "(if not set, assume first N atoms)"));
374 settings.impl_->plotSettings.initOptions(options);
376 // Add common options for trajectory processing.
377 if (!settings.hasFlag(TrajectoryAnalysisSettings::efNoUserRmPBC))
379 options->addOption(BooleanOption("rmpbc").store(&settings.impl_->bRmPBC)
380 .description("Make molecules whole for each frame"));
382 if (!settings.hasFlag(TrajectoryAnalysisSettings::efNoUserPBC))
384 options->addOption(BooleanOption("pbc").store(&settings.impl_->bPBC)
385 .description("Use periodic boundary conditions for distance calculation"));
391 TrajectoryAnalysisRunnerCommon::optionsFinished()
393 if (impl_->trjfile_.empty() && impl_->topfile_.empty())
395 GMX_THROW(InconsistentInputError("No trajectory or topology provided, nothing to do!"));
398 if (impl_->trajectoryGroup_.isValid() && impl_->trjfile_.empty())
400 GMX_THROW(InconsistentInputError("-fgroup only makes sense together with a trajectory (-f)"));
403 impl_->settings_.impl_->plotSettings.setTimeUnit(impl_->settings_.timeUnit());
405 if (impl_->bStartTimeSet_)
407 setTimeValue(TBEGIN, impl_->startTime_);
409 if (impl_->bEndTimeSet_)
411 setTimeValue(TEND, impl_->endTime_);
413 if (impl_->bDeltaTimeSet_)
415 setTimeValue(TDELTA, impl_->deltaTime_);
421 TrajectoryAnalysisRunnerCommon::initTopology()
423 const bool topologyRequired =
424 impl_->settings_.hasFlag(TrajectoryAnalysisSettings::efRequireTop);
425 impl_->initTopology(topologyRequired);
430 TrajectoryAnalysisRunnerCommon::initFirstFrame()
432 impl_->initFirstFrame();
437 TrajectoryAnalysisRunnerCommon::initFrameIndexGroup()
439 impl_->initFrameIndexGroup();
444 TrajectoryAnalysisRunnerCommon::readNextFrame()
446 bool bContinue = false;
449 bContinue = read_next_frame(impl_->oenv_, impl_->status_, impl_->fr);
453 impl_->finishTrajectory();
460 TrajectoryAnalysisRunnerCommon::initFrame()
462 if (impl_->gpbc_ != nullptr)
464 gmx_rmpbc_trxfr(impl_->gpbc_, impl_->fr);
470 TrajectoryAnalysisRunnerCommon::hasTrajectory() const
472 return impl_->hasTrajectory();
476 const TopologyInformation &
477 TrajectoryAnalysisRunnerCommon::topologyInformation() const
479 return impl_->topInfo_;
484 TrajectoryAnalysisRunnerCommon::frame() const
486 GMX_RELEASE_ASSERT(impl_->fr != nullptr, "Frame not available when accessed");