2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2010,2011,2012,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.
37 * Implements gmx::TrajectoryAnalysisRunnerCommon.
39 * \author Teemu Murtola <teemu.murtola@gmail.com>
40 * \ingroup module_trajectoryanalysis
44 #include "runnercommon.h"
51 #include "gromacs/compat/make_unique.h"
52 #include "gromacs/fileio/confio.h"
53 #include "gromacs/fileio/oenv.h"
54 #include "gromacs/fileio/timecontrol.h"
55 #include "gromacs/fileio/trxio.h"
56 #include "gromacs/math/vec.h"
57 #include "gromacs/options/basicoptions.h"
58 #include "gromacs/options/filenameoption.h"
59 #include "gromacs/options/ioptionscontainer.h"
60 #include "gromacs/pbcutil/rmpbc.h"
61 #include "gromacs/selection/selection.h"
62 #include "gromacs/selection/selectioncollection.h"
63 #include "gromacs/selection/selectionoption.h"
64 #include "gromacs/selection/selectionoptionbehavior.h"
65 #include "gromacs/topology/topology.h"
66 #include "gromacs/trajectory/trajectoryframe.h"
67 #include "gromacs/trajectoryanalysis/analysissettings.h"
68 #include "gromacs/utility/cstringutil.h"
69 #include "gromacs/utility/exceptions.h"
70 #include "gromacs/utility/gmxassert.h"
71 #include "gromacs/utility/programcontext.h"
72 #include "gromacs/utility/smalloc.h"
73 #include "gromacs/utility/stringutil.h"
75 #include "analysissettings-impl.h"
80 class TrajectoryAnalysisRunnerCommon::Impl
: public ITopologyProvider
83 explicit Impl(TrajectoryAnalysisSettings
*settings
);
86 bool hasTrajectory() const { return !trjfile_
.empty(); }
88 void initTopology(bool required
);
89 void initFirstFrame();
90 void initFrameIndexGroup();
91 void finishTrajectory();
93 // From ITopologyProvider
94 virtual gmx_mtop_t
*getTopology(bool required
)
96 initTopology(required
);
97 return topInfo_
.mtop_
.get();
99 virtual int getAtomCount()
101 if (!topInfo_
.hasTopology())
103 if (trajectoryGroup_
.isValid())
105 GMX_THROW(InconsistentInputError("-fgroup is only supported when -s is also specified"));
107 // Read the first frame if we don't know the maximum number of
115 TrajectoryAnalysisSettings
&settings_
;
116 TopologyInformation topInfo_
;
118 //! Name of the trajectory file (empty if not provided).
119 std::string trjfile_
;
120 //! Name of the topology file (empty if no topology provided).
121 std::string topfile_
;
122 Selection trajectoryGroup_
;
131 //! The current frame, or \p NULL if no frame loaded yet.
134 //! Used to store the status variable from read_first_frame().
135 t_trxstatus
*status_
;
136 gmx_output_env_t
*oenv_
;
140 TrajectoryAnalysisRunnerCommon::Impl::Impl(TrajectoryAnalysisSettings
*settings
)
141 : settings_(*settings
),
142 startTime_(0.0), endTime_(0.0), deltaTime_(0.0),
143 bStartTimeSet_(false), bEndTimeSet_(false), bDeltaTimeSet_(false),
144 bTrajOpen_(false), fr(nullptr), gpbc_(nullptr), status_(nullptr), oenv_(nullptr)
149 TrajectoryAnalysisRunnerCommon::Impl::~Impl()
154 // There doesn't seem to be a function for freeing frame data
161 if (oenv_
!= nullptr)
163 output_env_done(oenv_
);
168 TrajectoryAnalysisRunnerCommon::Impl::initTopology(bool required
)
170 // Return immediately if the topology has already been loaded.
171 if (topInfo_
.hasTopology())
176 if (required
&& topfile_
.empty())
178 GMX_THROW(InconsistentInputError("No topology provided, but one is required for analysis"));
181 // Load the topology if requested.
182 if (!topfile_
.empty())
184 topInfo_
.mtop_
= gmx::compat::make_unique
<gmx_mtop_t
>();
185 readConfAndTopology(topfile_
.c_str(), &topInfo_
.bTop_
, topInfo_
.mtop_
.get(),
186 &topInfo_
.ePBC_
, &topInfo_
.xtop_
, nullptr,
188 // TODO: Only load this here if the tool actually needs it; selections
189 // take care of themselves.
190 for (gmx_moltype_t
&moltype
: topInfo_
.mtop_
->moltype
)
192 if (!moltype
.atoms
.haveMass
)
194 // Try to read masses from database, be silent about missing masses
195 atomsSetMassesBasedOnNames(&moltype
.atoms
, FALSE
);
199 && !settings_
.hasFlag(TrajectoryAnalysisSettings::efUseTopX
))
201 sfree(topInfo_
.xtop_
);
202 topInfo_
.xtop_
= nullptr;
208 TrajectoryAnalysisRunnerCommon::Impl::initFirstFrame()
210 // Return if we have already initialized the trajectory.
215 time_unit_t time_unit
216 = static_cast<time_unit_t
>(settings_
.timeUnit() + 1);
217 output_env_init(&oenv_
, getProgramContext(), time_unit
, FALSE
, exvgNONE
, 0);
219 int frflags
= settings_
.frflags();
220 frflags
|= TRX_NEED_X
;
226 if (!read_first_frame(oenv_
, &status_
, trjfile_
.c_str(), fr
, frflags
))
228 GMX_THROW(FileIOError("Could not read coordinates from trajectory"));
232 if (topInfo_
.hasTopology())
234 const int topologyAtomCount
= topInfo_
.topology()->atoms
.nr
;
235 if (fr
->natoms
> topologyAtomCount
)
237 const std::string message
238 = formatString("Trajectory (%d atoms) does not match topology (%d atoms)",
239 fr
->natoms
, topologyAtomCount
);
240 GMX_THROW(InconsistentInputError(message
));
246 // Prepare a frame from topology information.
247 // TODO: Initialize more of the fields.
248 if (frflags
& (TRX_NEED_V
))
250 GMX_THROW(NotImplementedError("Velocity reading from a topology not implemented"));
252 if (frflags
& (TRX_NEED_F
))
254 GMX_THROW(InvalidInputError("Forces cannot be read from a topology"));
256 fr
->natoms
= topInfo_
.topology()->atoms
.nr
;
258 snew(fr
->x
, fr
->natoms
);
259 memcpy(fr
->x
, topInfo_
.xtop_
,
260 sizeof(*fr
->x
) * fr
->natoms
);
262 copy_mat(topInfo_
.boxtop_
, fr
->box
);
265 set_trxframe_ePBC(fr
, topInfo_
.ePBC());
266 if (topInfo_
.hasTopology() && settings_
.hasRmPBC())
268 gpbc_
= gmx_rmpbc_init(&topInfo_
.topology()->idef
, topInfo_
.ePBC(),
274 TrajectoryAnalysisRunnerCommon::Impl::initFrameIndexGroup()
276 if (!trajectoryGroup_
.isValid())
280 GMX_RELEASE_ASSERT(bTrajOpen_
,
281 "Trajectory index only makes sense with a real trajectory");
282 if (trajectoryGroup_
.atomCount() != fr
->natoms
)
284 const std::string message
= formatString(
285 "Selection specified with -fgroup has %d atoms, but "
286 "the trajectory (-f) has %d atoms.",
287 trajectoryGroup_
.atomCount(), fr
->natoms
);
288 GMX_THROW(InconsistentInputError(message
));
291 snew(fr
->index
, trajectoryGroup_
.atomCount());
292 std::copy(trajectoryGroup_
.atomIndices().begin(),
293 trajectoryGroup_
.atomIndices().end(),
298 TrajectoryAnalysisRunnerCommon::Impl::finishTrajectory()
305 if (gpbc_
!= nullptr)
307 gmx_rmpbc_done(gpbc_
);
312 /*********************************************************************
313 * TrajectoryAnalysisRunnerCommon
316 TrajectoryAnalysisRunnerCommon::TrajectoryAnalysisRunnerCommon(
317 TrajectoryAnalysisSettings
*settings
)
318 : impl_(new Impl(settings
))
323 TrajectoryAnalysisRunnerCommon::~TrajectoryAnalysisRunnerCommon()
329 TrajectoryAnalysisRunnerCommon::topologyProvider()
336 TrajectoryAnalysisRunnerCommon::initOptions(IOptionsContainer
*options
,
337 TimeUnitBehavior
*timeUnitBehavior
)
339 TrajectoryAnalysisSettings
&settings
= impl_
->settings_
;
341 // Add common file name arguments.
342 options
->addOption(FileNameOption("f")
343 .filetype(eftTrajectory
).inputFile()
344 .store(&impl_
->trjfile_
)
345 .defaultBasename("traj")
346 .description("Input trajectory or single configuration"));
347 options
->addOption(FileNameOption("s")
348 .filetype(eftTopology
).inputFile()
349 .store(&impl_
->topfile_
)
350 .defaultBasename("topol")
351 .description("Input structure"));
353 // Add options for trajectory time control.
354 options
->addOption(DoubleOption("b")
355 .store(&impl_
->startTime_
).storeIsSet(&impl_
->bStartTimeSet_
)
357 .description("First frame (%t) to read from trajectory"));
358 options
->addOption(DoubleOption("e")
359 .store(&impl_
->endTime_
).storeIsSet(&impl_
->bEndTimeSet_
)
361 .description("Last frame (%t) to read from trajectory"));
362 options
->addOption(DoubleOption("dt")
363 .store(&impl_
->deltaTime_
).storeIsSet(&impl_
->bDeltaTimeSet_
)
365 .description("Only use frame if t MOD dt == first time (%t)"));
367 // Add time unit option.
368 timeUnitBehavior
->setTimeUnitFromEnvironment();
369 timeUnitBehavior
->addTimeUnitOption(options
, "tu");
370 timeUnitBehavior
->setTimeUnitStore(&impl_
->settings_
.impl_
->timeUnit
);
372 options
->addOption(SelectionOption("fgroup")
373 .store(&impl_
->trajectoryGroup_
)
374 .onlySortedAtoms().onlyStatic()
375 .description("Atoms stored in the trajectory file "
376 "(if not set, assume first N atoms)"));
379 settings
.impl_
->plotSettings
.initOptions(options
);
381 // Add common options for trajectory processing.
382 if (!settings
.hasFlag(TrajectoryAnalysisSettings::efNoUserRmPBC
))
384 options
->addOption(BooleanOption("rmpbc").store(&settings
.impl_
->bRmPBC
)
385 .description("Make molecules whole for each frame"));
387 if (!settings
.hasFlag(TrajectoryAnalysisSettings::efNoUserPBC
))
389 options
->addOption(BooleanOption("pbc").store(&settings
.impl_
->bPBC
)
390 .description("Use periodic boundary conditions for distance calculation"));
396 TrajectoryAnalysisRunnerCommon::optionsFinished()
398 if (impl_
->trjfile_
.empty() && impl_
->topfile_
.empty())
400 GMX_THROW(InconsistentInputError("No trajectory or topology provided, nothing to do!"));
403 if (impl_
->trajectoryGroup_
.isValid() && impl_
->trjfile_
.empty())
405 GMX_THROW(InconsistentInputError("-fgroup only makes sense together with a trajectory (-f)"));
408 impl_
->settings_
.impl_
->plotSettings
.setTimeUnit(impl_
->settings_
.timeUnit());
410 if (impl_
->bStartTimeSet_
)
412 setTimeValue(TBEGIN
, impl_
->startTime_
);
414 if (impl_
->bEndTimeSet_
)
416 setTimeValue(TEND
, impl_
->endTime_
);
418 if (impl_
->bDeltaTimeSet_
)
420 setTimeValue(TDELTA
, impl_
->deltaTime_
);
426 TrajectoryAnalysisRunnerCommon::initTopology()
428 const bool topologyRequired
=
429 impl_
->settings_
.hasFlag(TrajectoryAnalysisSettings::efRequireTop
);
430 impl_
->initTopology(topologyRequired
);
435 TrajectoryAnalysisRunnerCommon::initFirstFrame()
437 impl_
->initFirstFrame();
442 TrajectoryAnalysisRunnerCommon::initFrameIndexGroup()
444 impl_
->initFrameIndexGroup();
449 TrajectoryAnalysisRunnerCommon::readNextFrame()
451 bool bContinue
= false;
454 bContinue
= read_next_frame(impl_
->oenv_
, impl_
->status_
, impl_
->fr
);
458 impl_
->finishTrajectory();
465 TrajectoryAnalysisRunnerCommon::initFrame()
467 if (impl_
->gpbc_
!= nullptr)
469 gmx_rmpbc_trxfr(impl_
->gpbc_
, impl_
->fr
);
475 TrajectoryAnalysisRunnerCommon::hasTrajectory() const
477 return impl_
->hasTrajectory();
481 const TopologyInformation
&
482 TrajectoryAnalysisRunnerCommon::topologyInformation() const
484 return impl_
->topInfo_
;
489 TrajectoryAnalysisRunnerCommon::frame() const
491 GMX_RELEASE_ASSERT(impl_
->fr
!= nullptr, "Frame not available when accessed");