2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2010,2011,2012,2013,2014,2015, 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"
48 #include "gromacs/fileio/confio.h"
49 #include "gromacs/fileio/oenv.h"
50 #include "gromacs/fileio/timecontrol.h"
51 #include "gromacs/fileio/trx.h"
52 #include "gromacs/fileio/trxio.h"
53 #include "gromacs/math/vec.h"
54 #include "gromacs/options/basicoptions.h"
55 #include "gromacs/options/filenameoption.h"
56 #include "gromacs/options/ioptionscontainer.h"
57 #include "gromacs/options/options.h"
58 #include "gromacs/pbcutil/rmpbc.h"
59 #include "gromacs/selection/indexutil.h"
60 #include "gromacs/selection/selectioncollection.h"
61 #include "gromacs/selection/selectionfileoption.h"
62 #include "gromacs/topology/topology.h"
63 #include "gromacs/trajectoryanalysis/analysissettings.h"
64 #include "gromacs/utility/cstringutil.h"
65 #include "gromacs/utility/exceptions.h"
66 #include "gromacs/utility/gmxassert.h"
67 #include "gromacs/utility/programcontext.h"
68 #include "gromacs/utility/smalloc.h"
69 #include "gromacs/utility/stringutil.h"
71 #include "analysissettings-impl.h"
76 class TrajectoryAnalysisRunnerCommon::Impl
79 Impl(TrajectoryAnalysisSettings
*settings
);
82 void finishTrajectory();
84 TrajectoryAnalysisSettings
&settings_
;
85 TopologyInformation topInfo_
;
87 //! Name of the trajectory file (empty if not provided).
89 //! Name of the topology file (empty if no topology provided).
91 //! Name of the index file (empty if no index file provided).
100 gmx_ana_indexgrps_t
*grps_
;
102 //! The current frame, or \p NULL if no frame loaded yet.
105 //! Used to store the status variable from read_first_frame().
106 t_trxstatus
*status_
;
107 gmx_output_env_t
*oenv_
;
111 TrajectoryAnalysisRunnerCommon::Impl::Impl(TrajectoryAnalysisSettings
*settings
)
112 : settings_(*settings
),
113 startTime_(0.0), endTime_(0.0), deltaTime_(0.0),
114 bStartTimeSet_(false), bEndTimeSet_(false), bDeltaTimeSet_(false),
116 bTrajOpen_(false), fr(NULL
), gpbc_(NULL
), status_(NULL
), oenv_(NULL
)
121 TrajectoryAnalysisRunnerCommon::Impl::~Impl()
125 gmx_ana_indexgrps_free(grps_
);
130 // There doesn't seem to be a function for freeing frame data
138 output_env_done(oenv_
);
144 TrajectoryAnalysisRunnerCommon::Impl::finishTrajectory()
153 gmx_rmpbc_done(gpbc_
);
158 /*********************************************************************
159 * TrajectoryAnalysisRunnerCommon
162 TrajectoryAnalysisRunnerCommon::TrajectoryAnalysisRunnerCommon(
163 TrajectoryAnalysisSettings
*settings
)
164 : impl_(new Impl(settings
))
169 TrajectoryAnalysisRunnerCommon::~TrajectoryAnalysisRunnerCommon()
175 TrajectoryAnalysisRunnerCommon::initOptions(IOptionsContainer
*options
,
176 TimeUnitBehavior
*timeUnitBehavior
)
178 TrajectoryAnalysisSettings
&settings
= impl_
->settings_
;
180 // Add common file name arguments.
181 options
->addOption(FileNameOption("f")
182 .filetype(eftTrajectory
).inputFile()
183 .store(&impl_
->trjfile_
)
184 .defaultBasename("traj")
185 .description("Input trajectory or single configuration"));
186 options
->addOption(FileNameOption("s")
187 .filetype(eftTopology
).inputFile()
188 .store(&impl_
->topfile_
)
189 .defaultBasename("topol")
190 .description("Input structure"));
191 options
->addOption(FileNameOption("n")
192 .filetype(eftIndex
).inputFile()
193 .store(&impl_
->ndxfile_
)
194 .defaultBasename("index")
195 .description("Extra index groups"));
197 // Add options for trajectory time control.
198 options
->addOption(DoubleOption("b")
199 .store(&impl_
->startTime_
).storeIsSet(&impl_
->bStartTimeSet_
)
201 .description("First frame (%t) to read from trajectory"));
202 options
->addOption(DoubleOption("e")
203 .store(&impl_
->endTime_
).storeIsSet(&impl_
->bEndTimeSet_
)
205 .description("Last frame (%t) to read from trajectory"));
206 options
->addOption(DoubleOption("dt")
207 .store(&impl_
->deltaTime_
).storeIsSet(&impl_
->bDeltaTimeSet_
)
209 .description("Only use frame if t MOD dt == first time (%t)"));
211 // Add time unit option.
212 timeUnitBehavior
->setTimeUnitFromEnvironment();
213 timeUnitBehavior
->addTimeUnitOption(options
, "tu");
214 timeUnitBehavior
->setTimeUnitStore(&impl_
->settings_
.impl_
->timeUnit
);
217 settings
.impl_
->plotSettings
.initOptions(options
);
219 // Add common options for trajectory processing.
220 if (!settings
.hasFlag(TrajectoryAnalysisSettings::efNoUserRmPBC
))
222 options
->addOption(BooleanOption("rmpbc").store(&settings
.impl_
->bRmPBC
)
223 .description("Make molecules whole for each frame"));
225 if (!settings
.hasFlag(TrajectoryAnalysisSettings::efNoUserPBC
))
227 options
->addOption(BooleanOption("pbc").store(&settings
.impl_
->bPBC
)
228 .description("Use periodic boundary conditions for distance calculation"));
231 options
->addOption(SelectionFileOption("sf"));
236 TrajectoryAnalysisRunnerCommon::optionsFinished()
238 impl_
->settings_
.impl_
->plotSettings
.setTimeUnit(impl_
->settings_
.timeUnit());
240 if (impl_
->trjfile_
.empty() && impl_
->topfile_
.empty())
242 GMX_THROW(InconsistentInputError("No trajectory or topology provided, nothing to do!"));
245 if (impl_
->bStartTimeSet_
)
247 setTimeValue(TBEGIN
, impl_
->startTime_
);
249 if (impl_
->bEndTimeSet_
)
251 setTimeValue(TEND
, impl_
->endTime_
);
253 if (impl_
->bDeltaTimeSet_
)
255 setTimeValue(TDELTA
, impl_
->deltaTime_
);
261 TrajectoryAnalysisRunnerCommon::initIndexGroups(SelectionCollection
*selections
,
264 if (impl_
->ndxfile_
.empty())
268 selections
->setIndexGroups(NULL
);
271 initTopology(selections
);
273 const char *const ndxfile
274 = (!impl_
->ndxfile_
.empty() ? impl_
->ndxfile_
.c_str() : NULL
);
275 gmx_ana_indexgrps_init(&impl_
->grps_
, impl_
->topInfo_
.topology(), ndxfile
);
276 selections
->setIndexGroups(impl_
->grps_
);
281 TrajectoryAnalysisRunnerCommon::doneIndexGroups(SelectionCollection
*selections
)
283 if (impl_
->grps_
!= NULL
)
285 selections
->setIndexGroups(NULL
);
286 gmx_ana_indexgrps_free(impl_
->grps_
);
293 TrajectoryAnalysisRunnerCommon::initTopology(SelectionCollection
*selections
)
295 // Return immediately if the topology has already been loaded.
296 if (impl_
->topInfo_
.hasTopology())
301 const TrajectoryAnalysisSettings
&settings
= impl_
->settings_
;
302 const bool bRequireTop
303 = settings
.hasFlag(TrajectoryAnalysisSettings::efRequireTop
)
304 || selections
->requiresTopology();
305 if (bRequireTop
&& impl_
->topfile_
.empty())
307 GMX_THROW(InconsistentInputError("No topology provided, but one is required for analysis"));
310 // Load the topology if requested.
311 if (!impl_
->topfile_
.empty())
313 snew(impl_
->topInfo_
.top_
, 1);
314 impl_
->topInfo_
.bTop_
= read_tps_conf(impl_
->topfile_
.c_str(),
315 impl_
->topInfo_
.top_
, &impl_
->topInfo_
.ePBC_
,
316 &impl_
->topInfo_
.xtop_
, NULL
, impl_
->topInfo_
.boxtop_
, TRUE
);
318 && !settings
.hasFlag(TrajectoryAnalysisSettings::efUseTopX
))
320 sfree(impl_
->topInfo_
.xtop_
);
321 impl_
->topInfo_
.xtop_
= NULL
;
325 // Read the first frame if we don't know the maximum number of atoms
328 if (!impl_
->topInfo_
.hasTopology())
331 natoms
= impl_
->fr
->natoms
;
333 selections
->setTopology(impl_
->topInfo_
.topology(), natoms
);
338 gmx_ana_poscalc_coll_print_tree(stderr, impl_->pcc);
339 fprintf(stderr, "\n");
346 TrajectoryAnalysisRunnerCommon::initFirstFrame()
348 // Return if we have already initialized the trajectory.
353 time_unit_t time_unit
354 = static_cast<time_unit_t
>(impl_
->settings_
.timeUnit() + 1);
355 output_env_init(&impl_
->oenv_
, getProgramContext(), time_unit
, FALSE
, exvgNONE
, 0);
357 int frflags
= impl_
->settings_
.frflags();
358 frflags
|= TRX_NEED_X
;
362 const TopologyInformation
&top
= impl_
->topInfo_
;
365 if (!read_first_frame(impl_
->oenv_
, &impl_
->status_
,
366 impl_
->trjfile_
.c_str(), impl_
->fr
, frflags
))
368 GMX_THROW(FileIOError("Could not read coordinates from trajectory"));
370 impl_
->bTrajOpen_
= true;
372 if (top
.hasTopology() && impl_
->fr
->natoms
> top
.topology()->atoms
.nr
)
374 GMX_THROW(InconsistentInputError(formatString(
375 "Trajectory (%d atoms) does not match topology (%d atoms)",
376 impl_
->fr
->natoms
, top
.topology()->atoms
.nr
)));
381 // Prepare a frame from topology information.
382 // TODO: Initialize more of the fields.
383 if (frflags
& (TRX_NEED_V
))
385 GMX_THROW(NotImplementedError("Velocity reading from a topology not implemented"));
387 if (frflags
& (TRX_NEED_F
))
389 GMX_THROW(InvalidInputError("Forces cannot be read from a topology"));
391 impl_
->fr
->flags
= frflags
;
392 impl_
->fr
->natoms
= top
.topology()->atoms
.nr
;
393 impl_
->fr
->bX
= TRUE
;
394 snew(impl_
->fr
->x
, impl_
->fr
->natoms
);
395 memcpy(impl_
->fr
->x
, top
.xtop_
,
396 sizeof(*impl_
->fr
->x
) * impl_
->fr
->natoms
);
397 impl_
->fr
->bBox
= TRUE
;
398 copy_mat(const_cast<rvec
*>(top
.boxtop_
), impl_
->fr
->box
);
401 set_trxframe_ePBC(impl_
->fr
, top
.ePBC());
402 if (top
.hasTopology() && impl_
->settings_
.hasRmPBC())
404 impl_
->gpbc_
= gmx_rmpbc_init(&top
.topology()->idef
, top
.ePBC(),
411 TrajectoryAnalysisRunnerCommon::readNextFrame()
413 bool bContinue
= false;
416 bContinue
= read_next_frame(impl_
->oenv_
, impl_
->status_
, impl_
->fr
);
420 impl_
->finishTrajectory();
427 TrajectoryAnalysisRunnerCommon::initFrame()
429 if (impl_
->gpbc_
!= NULL
)
431 gmx_rmpbc_trxfr(impl_
->gpbc_
, impl_
->fr
);
437 TrajectoryAnalysisRunnerCommon::hasTrajectory() const
439 return !impl_
->trjfile_
.empty();
443 const TopologyInformation
&
444 TrajectoryAnalysisRunnerCommon::topologyInformation() const
446 return impl_
->topInfo_
;
451 TrajectoryAnalysisRunnerCommon::frame() const
453 GMX_RELEASE_ASSERT(impl_
->fr
!= NULL
, "Frame not available when accessed");