2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2016,2018,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.
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::analysismodules::Trajectory.
39 * \author Teemu Murtola <teemu.murtola@gmail.com>
40 * \ingroup module_trajectoryanalysis
44 #include "trajectory.h"
48 #include "gromacs/analysisdata/analysisdata.h"
49 #include "gromacs/analysisdata/modules/plot.h"
50 #include "gromacs/fileio/trxio.h"
51 #include "gromacs/options/basicoptions.h"
52 #include "gromacs/options/filenameoption.h"
53 #include "gromacs/options/ioptionscontainer.h"
54 #include "gromacs/selection/selection.h"
55 #include "gromacs/selection/selectionoption.h"
56 #include "gromacs/trajectory/trajectoryframe.h"
57 #include "gromacs/trajectoryanalysis/analysissettings.h"
62 namespace analysismodules
68 /********************************************************************
72 class Trajectory : public TrajectoryAnalysisModule
77 void initOptions(IOptionsContainer* options, TrajectoryAnalysisSettings* settings) override;
78 void optionsFinished(TrajectoryAnalysisSettings* settings) override;
79 void initAnalysis(const TrajectoryAnalysisSettings& settings, const TopologyInformation& top) override;
81 void analyzeFrame(int frnr, const t_trxframe& fr, t_pbc* pbc, TrajectoryAnalysisModuleData* pdata) override;
83 void finishAnalysis(int nframes) override;
84 void writeOutput() override;
92 std::array<bool, 4> dimMask_;
93 std::array<bool, 4> maskSet_;
100 Trajectory::Trajectory() : dimMask_{ true, true, true, false }, maskSet_{}
102 registerAnalysisDataset(&xdata_, "x");
103 registerAnalysisDataset(&vdata_, "v");
104 registerAnalysisDataset(&fdata_, "f");
108 void Trajectory::initOptions(IOptionsContainer* options, TrajectoryAnalysisSettings* settings)
110 static const char* const desc[] = {
111 "[THISMODULE] plots coordinates, velocities, and/or forces for",
112 "provided selections. By default, the X, Y, and Z components for",
113 "the requested vectors are plotted, but specifying one or more of",
114 "[TT]-len[tt], [TT]-x[tt], [TT]-y[tt], and [TT]-z[tt] overrides this.",
116 "For dynamic selections, currently the values are written out for",
117 "all positions that the selection could select."
120 settings->setHelpText(desc);
122 options->addOption(FileNameOption("ox")
123 .filetype(OptionFileType::Plot)
126 .defaultBasename("coord")
127 .description("Coordinates for each position as a function of time"));
128 options->addOption(FileNameOption("ov")
129 .filetype(OptionFileType::Plot)
132 .defaultBasename("veloc")
133 .description("Velocities for each position as a function of time"));
134 options->addOption(FileNameOption("of")
135 .filetype(OptionFileType::Plot)
138 .defaultBasename("force")
139 .description("Forces for each position as a function of time"));
142 SelectionOption("select").storeVector(&sel_).required().dynamicMask().multiValue().description(
143 "Selections to analyze"));
146 BooleanOption("x").store(&dimMask_[XX]).storeIsSet(&maskSet_[XX]).description("Plot X component"));
148 BooleanOption("y").store(&dimMask_[YY]).storeIsSet(&maskSet_[YY]).description("Plot Y component"));
150 BooleanOption("z").store(&dimMask_[ZZ]).storeIsSet(&maskSet_[ZZ]).description("Plot Z component"));
152 BooleanOption("len").store(&dimMask_[DIM]).storeIsSet(&maskSet_[DIM]).description("Plot vector length"));
156 void Trajectory::optionsFinished(TrajectoryAnalysisSettings* settings)
158 int frameFlags = TRX_NEED_X;
161 frameFlags |= TRX_READ_V;
165 frameFlags |= TRX_READ_F;
167 settings->setFrameFlags(frameFlags);
168 if (std::count(std::begin(maskSet_), std::end(maskSet_), true) > 0)
170 for (int i = 0; i <= DIM; ++i)
181 void Trajectory::initAnalysis(const TrajectoryAnalysisSettings& settings, const TopologyInformation& /*top*/)
185 xdata_.setDataSetCount(sel_.size());
186 for (size_t g = 0; g < sel_.size(); ++g)
188 xdata_.setColumnCount(g, 3 * sel_[g].posCount());
190 AnalysisDataVectorPlotModulePointer plot(new AnalysisDataVectorPlotModule(settings.plotSettings()));
191 plot->setWriteMask(dimMask_.data());
192 plot->setFileName(fnX_);
193 plot->setTitle("Coordinates");
194 plot->setXAxisIsTime();
195 plot->setYLabel("Value [nm]");
196 xdata_.addModule(plot);
200 vdata_.setDataSetCount(sel_.size());
201 for (size_t g = 0; g < sel_.size(); ++g)
203 sel_[g].setEvaluateVelocities(true);
204 vdata_.setColumnCount(g, 3 * sel_[g].posCount());
206 AnalysisDataVectorPlotModulePointer plot(new AnalysisDataVectorPlotModule(settings.plotSettings()));
207 plot->setWriteMask(dimMask_.data());
208 plot->setFileName(fnV_);
209 plot->setTitle("Velocities");
210 plot->setXAxisIsTime();
211 plot->setYLabel("Value [nm/ps]");
212 vdata_.addModule(plot);
216 fdata_.setDataSetCount(sel_.size());
217 for (size_t g = 0; g < sel_.size(); ++g)
219 sel_[g].setEvaluateForces(true);
220 fdata_.setColumnCount(g, 3 * sel_[g].posCount());
222 AnalysisDataVectorPlotModulePointer plot(new AnalysisDataVectorPlotModule(settings.plotSettings()));
223 plot->setWriteMask(dimMask_.data());
224 plot->setFileName(fnF_);
225 plot->setTitle("Forces");
226 plot->setXAxisIsTime();
227 plot->setYLabel("Value [kJ mol\\S-1\\N nm\\S-1\\N]");
228 fdata_.addModule(plot);
233 //! Helper function for Trajectory::analyzeFrame
235 void analyzeFrameImpl(int frnr, const t_trxframe& fr, AnalysisDataHandle* dh, const SelectionList& sel, T getField)
239 dh->startFrame(frnr, fr.time);
240 for (size_t g = 0; g < sel.size(); ++g)
242 dh->selectDataSet(g);
243 for (int i = 0; i < sel[g].posCount(); ++i)
245 const SelectionPosition& pos = sel[g].position(i);
246 dh->setPoints(i * 3, 3, getField(pos), pos.selected());
253 void Trajectory::analyzeFrame(int frnr, const t_trxframe& fr, t_pbc* /* pbc */, TrajectoryAnalysisModuleData* pdata)
255 AnalysisDataHandle dh = pdata->dataHandle(xdata_);
256 const SelectionList& sel = TrajectoryAnalysisModuleData::parallelSelections(sel_);
257 analyzeFrameImpl(frnr, fr, &dh, sel, [](const SelectionPosition& pos) { return pos.x(); });
260 analyzeFrameImpl(frnr, fr, &dh, sel, [](const SelectionPosition& pos) { return pos.v(); });
264 analyzeFrameImpl(frnr, fr, &dh, sel, [](const SelectionPosition& pos) { return pos.f(); });
269 void Trajectory::finishAnalysis(int /*nframes*/) {}
272 void Trajectory::writeOutput() {}
276 const char TrajectoryInfo::name[] = "trajectory";
277 const char TrajectoryInfo::shortDescription[] =
278 "Print coordinates, velocities, and/or forces for selections";
280 TrajectoryAnalysisModulePointer TrajectoryInfo::create()
282 return TrajectoryAnalysisModulePointer(new Trajectory);
285 } // namespace analysismodules