2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2016,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::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 virtual void initOptions(IOptionsContainer *options,
78 TrajectoryAnalysisSettings *settings);
79 virtual void optionsFinished(TrajectoryAnalysisSettings *settings);
80 virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
81 const TopologyInformation &top);
83 virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
84 TrajectoryAnalysisModuleData *pdata);
86 virtual void finishAnalysis(int nframes);
87 virtual void writeOutput();
95 std::array<bool, 4> dimMask_;
96 std::array<bool, 4> maskSet_;
103 Trajectory::Trajectory() :
104 dimMask_ {true, true, true, false}, maskSet_ {}
106 registerAnalysisDataset(&xdata_, "x");
107 registerAnalysisDataset(&vdata_, "v");
108 registerAnalysisDataset(&fdata_, "f");
113 Trajectory::initOptions(IOptionsContainer *options, TrajectoryAnalysisSettings *settings)
115 static const char *const desc[] = {
116 "[THISMODULE] plots coordinates, velocities, and/or forces for",
117 "provided selections. By default, the X, Y, and Z components for",
118 "the requested vectors are plotted, but specifying one or more of",
119 "[TT]-len[tt], [TT]-x[tt], [TT]-y[tt], and [TT]-z[tt] overrides this.",
121 "For dynamic selections, currently the values are written out for",
122 "all positions that the selection could select."
125 settings->setHelpText(desc);
127 options->addOption(FileNameOption("ox").filetype(eftPlot).outputFile()
128 .store(&fnX_).defaultBasename("coord")
129 .description("Coordinates for each position as a function of time"));
130 options->addOption(FileNameOption("ov").filetype(eftPlot).outputFile()
131 .store(&fnV_).defaultBasename("veloc")
132 .description("Velocities for each position as a function of time"));
133 options->addOption(FileNameOption("of").filetype(eftPlot).outputFile()
134 .store(&fnF_).defaultBasename("force")
135 .description("Forces for each position as a function of time"));
137 options->addOption(SelectionOption("select").storeVector(&sel_)
138 .required().dynamicMask().multiValue()
139 .description("Selections to analyze"));
141 options->addOption(BooleanOption("x").store(&dimMask_[XX])
142 .storeIsSet(&maskSet_[XX])
143 .description("Plot X component"));
144 options->addOption(BooleanOption("y").store(&dimMask_[YY])
145 .storeIsSet(&maskSet_[YY])
146 .description("Plot Y component"));
147 options->addOption(BooleanOption("z").store(&dimMask_[ZZ])
148 .storeIsSet(&maskSet_[ZZ])
149 .description("Plot Z component"));
150 options->addOption(BooleanOption("len").store(&dimMask_[DIM])
151 .storeIsSet(&maskSet_[DIM])
152 .description("Plot vector length"));
157 Trajectory::optionsFinished(TrajectoryAnalysisSettings *settings)
159 int frameFlags = TRX_NEED_X;
162 frameFlags |= TRX_READ_V;
166 frameFlags |= TRX_READ_F;
168 settings->setFrameFlags(frameFlags);
169 if (std::count(std::begin(maskSet_), std::end(maskSet_), true) > 0)
171 for (int i = 0; i <= DIM; ++i)
183 Trajectory::initAnalysis(const TrajectoryAnalysisSettings &settings,
184 const TopologyInformation & /*top*/)
188 xdata_.setDataSetCount(sel_.size());
189 for (size_t g = 0; g < sel_.size(); ++g)
191 xdata_.setColumnCount(g, 3*sel_[g].posCount());
193 AnalysisDataVectorPlotModulePointer plot(
194 new AnalysisDataVectorPlotModule(settings.plotSettings()));
195 plot->setWriteMask(dimMask_.data());
196 plot->setFileName(fnX_);
197 plot->setTitle("Coordinates");
198 plot->setXAxisIsTime();
199 plot->setYLabel("Value [nm]");
200 xdata_.addModule(plot);
204 vdata_.setDataSetCount(sel_.size());
205 for (size_t g = 0; g < sel_.size(); ++g)
207 sel_[g].setEvaluateVelocities(true);
208 vdata_.setColumnCount(g, 3*sel_[g].posCount());
210 AnalysisDataVectorPlotModulePointer plot(
211 new AnalysisDataVectorPlotModule(settings.plotSettings()));
212 plot->setWriteMask(dimMask_.data());
213 plot->setFileName(fnV_);
214 plot->setTitle("Velocities");
215 plot->setXAxisIsTime();
216 plot->setYLabel("Value [nm/ps]");
217 vdata_.addModule(plot);
221 fdata_.setDataSetCount(sel_.size());
222 for (size_t g = 0; g < sel_.size(); ++g)
224 sel_[g].setEvaluateForces(true);
225 fdata_.setColumnCount(g, 3*sel_[g].posCount());
227 AnalysisDataVectorPlotModulePointer plot(
228 new AnalysisDataVectorPlotModule(settings.plotSettings()));
229 plot->setWriteMask(dimMask_.data());
230 plot->setFileName(fnF_);
231 plot->setTitle("Forces");
232 plot->setXAxisIsTime();
233 plot->setYLabel("Value [kJ mol\\S-1\\N nm\\S-1\\N]");
234 fdata_.addModule(plot);
239 //! Helper function for Trajectory::analyzeFrame
240 template<typename T> void
241 analyzeFrameImpl(int frnr, const t_trxframe &fr,
242 AnalysisDataHandle* dh, const SelectionList &sel, T getField)
246 dh->startFrame(frnr, fr.time);
247 for (size_t g = 0; g < sel.size(); ++g)
249 dh->selectDataSet(g);
250 for (int i = 0; i < sel[g].posCount(); ++i)
252 const SelectionPosition &pos = sel[g].position(i);
253 dh->setPoints(i*3, 3, getField(pos), pos.selected());
261 Trajectory::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc * /* pbc */,
262 TrajectoryAnalysisModuleData *pdata)
264 AnalysisDataHandle dh = pdata->dataHandle(xdata_);
265 const SelectionList &sel = pdata->parallelSelections(sel_);
266 analyzeFrameImpl(frnr, fr, &dh, sel, [](const SelectionPosition &pos) { return pos.x(); });
269 analyzeFrameImpl(frnr, fr, &dh, sel, [](const SelectionPosition &pos) { return pos.v(); });
273 analyzeFrameImpl(frnr, fr, &dh, sel, [](const SelectionPosition &pos) { return pos.f(); });
280 Trajectory::finishAnalysis(int /*nframes*/)
286 Trajectory::writeOutput()
292 const char TrajectoryInfo::name[] = "trajectory";
293 const char TrajectoryInfo::shortDescription[] =
294 "Print coordinates, velocities, and/or forces for selections";
296 TrajectoryAnalysisModulePointer TrajectoryInfo::create()
298 return TrajectoryAnalysisModulePointer(new Trajectory);
301 } // namespace analysismodules