2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2010,2011,2012,2013,2014, 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::Distance.
39 * \author Teemu Murtola <teemu.murtola@gmail.com>
40 * \ingroup module_trajectoryanalysis
46 #include "gromacs/analysisdata/analysisdata.h"
47 #include "gromacs/analysisdata/modules/average.h"
48 #include "gromacs/analysisdata/modules/histogram.h"
49 #include "gromacs/analysisdata/modules/plot.h"
50 #include "gromacs/fileio/trx.h"
51 #include "gromacs/math/vec.h"
52 #include "gromacs/options/basicoptions.h"
53 #include "gromacs/options/filenameoption.h"
54 #include "gromacs/options/options.h"
55 #include "gromacs/pbcutil/pbc.h"
56 #include "gromacs/selection/selection.h"
57 #include "gromacs/selection/selectionoption.h"
58 #include "gromacs/trajectoryanalysis/analysissettings.h"
59 #include "gromacs/utility/arrayref.h"
60 #include "gromacs/utility/exceptions.h"
61 #include "gromacs/utility/stringutil.h"
66 namespace analysismodules
72 class Distance
: public TrajectoryAnalysisModule
77 virtual void initOptions(Options
*options
,
78 TrajectoryAnalysisSettings
*settings
);
79 virtual void initAnalysis(const TrajectoryAnalysisSettings
&settings
,
80 const TopologyInformation
&top
);
82 virtual void analyzeFrame(int frnr
, const t_trxframe
&fr
, t_pbc
*pbc
,
83 TrajectoryAnalysisModuleData
*pdata
);
85 virtual void finishAnalysis(int nframes
);
86 virtual void writeOutput();
90 std::string fnAverage_
;
93 std::string fnHistogram_
;
94 std::string fnAllStats_
;
99 AnalysisData distances_
;
101 AnalysisDataAverageModulePointer summaryStatsModule_
;
102 AnalysisDataAverageModulePointer allStatsModule_
;
103 AnalysisDataFrameAverageModulePointer averageModule_
;
104 AnalysisDataSimpleHistogramModulePointer histogramModule_
;
106 // Copy and assign disallowed by base.
110 : TrajectoryAnalysisModule(DistanceInfo::name
, DistanceInfo::shortDescription
),
111 meanLength_(0.1), lengthDev_(1.0), binWidth_(0.001)
113 summaryStatsModule_
.reset(new AnalysisDataAverageModule());
114 summaryStatsModule_
->setAverageDataSets(true);
115 distances_
.addModule(summaryStatsModule_
);
116 allStatsModule_
.reset(new AnalysisDataAverageModule());
117 distances_
.addModule(allStatsModule_
);
118 averageModule_
.reset(new AnalysisDataFrameAverageModule());
119 distances_
.addModule(averageModule_
);
120 histogramModule_
.reset(new AnalysisDataSimpleHistogramModule());
121 distances_
.addModule(histogramModule_
);
123 registerAnalysisDataset(&distances_
, "dist");
124 registerAnalysisDataset(&xyz_
, "xyz");
125 registerBasicDataset(summaryStatsModule_
.get(), "stats");
126 registerBasicDataset(allStatsModule_
.get(), "allstats");
127 registerBasicDataset(averageModule_
.get(), "average");
128 registerBasicDataset(&histogramModule_
->averager(), "histogram");
133 Distance::initOptions(Options
*options
, TrajectoryAnalysisSettings
* /*settings*/)
135 static const char *const desc
[] = {
136 "[THISMODULE] calculates distances between pairs of positions",
137 "as a function of time. Each selection specifies an independent set",
138 "of distances to calculate. Each selection should consist of pairs",
139 "of positions, and the distances are computed between positions 1-2,",
141 "[TT]-oav[tt] writes the average distance as a function of time for",
143 "[TT]-oall[tt] writes all the individual distances.",
144 "[TT]-oxyz[tt] does the same, but the x, y, and z components of the",
145 "distance are written instead of the norm.",
146 "[TT]-oh[tt] writes a histogram of the distances for each selection.",
147 "The location of the histogram is set with [TT]-len[tt] and",
148 "[TT]-tol[tt]. Bin width is set with [TT]-binw[tt].",
149 "[TT]-oallstat[tt] writes out the average and standard deviation for",
150 "each individual distance, calculated over the frames."
153 options
->setDescription(desc
);
155 options
->addOption(FileNameOption("oav").filetype(eftPlot
).outputFile()
156 .store(&fnAverage_
).defaultBasename("distave")
157 .description("Average distances as function of time"));
158 options
->addOption(FileNameOption("oall").filetype(eftPlot
).outputFile()
159 .store(&fnAll_
).defaultBasename("dist")
160 .description("All distances as function of time"));
161 options
->addOption(FileNameOption("oxyz").filetype(eftPlot
).outputFile()
162 .store(&fnXYZ_
).defaultBasename("distxyz")
163 .description("Distance components as function of time"));
164 options
->addOption(FileNameOption("oh").filetype(eftPlot
).outputFile()
165 .store(&fnHistogram_
).defaultBasename("disthist")
166 .description("Histogram of the distances"));
167 options
->addOption(FileNameOption("oallstat").filetype(eftPlot
).outputFile()
168 .store(&fnAllStats_
).defaultBasename("diststat")
169 .description("Statistics for individual distances"));
170 options
->addOption(SelectionOption("select").storeVector(&sel_
)
171 .required().dynamicMask().multiValue()
172 .description("Position pairs to calculate distances for"));
173 // TODO: Extend the histogramming implementation to allow automatic
174 // extension of the histograms to cover the data, removing the need for
175 // the first two options.
176 options
->addOption(DoubleOption("len").store(&meanLength_
)
177 .description("Mean distance for histogramming"));
178 options
->addOption(DoubleOption("tol").store(&lengthDev_
)
179 .description("Width of full distribution as fraction of [TT]-len[tt]"));
180 options
->addOption(DoubleOption("binw").store(&binWidth_
)
181 .description("Bin width for histogramming"));
186 * Checks that selections conform to the expectations of the tool.
188 void checkSelections(const SelectionList
&sel
)
190 for (size_t g
= 0; g
< sel
.size(); ++g
)
192 if (sel
[g
].posCount() % 2 != 0)
194 std::string message
= formatString(
195 "Selection '%s' does not evaluate into an even number of positions "
196 "(there are %d positions)",
197 sel
[g
].name(), sel
[g
].posCount());
198 GMX_THROW(InconsistentInputError(message
));
200 if (sel
[g
].isDynamic())
202 for (int i
= 0; i
< sel
[g
].posCount(); i
+= 2)
204 if (sel
[g
].position(i
).selected() != sel
[g
].position(i
+1).selected())
206 std::string message
=
207 formatString("Dynamic selection %d does not select "
208 "a consistent set of pairs over the frames",
209 static_cast<int>(g
+ 1));
210 GMX_THROW(InconsistentInputError(message
));
219 Distance::initAnalysis(const TrajectoryAnalysisSettings
&settings
,
220 const TopologyInformation
& /*top*/)
222 checkSelections(sel_
);
224 distances_
.setDataSetCount(sel_
.size());
225 xyz_
.setDataSetCount(sel_
.size());
226 for (size_t i
= 0; i
< sel_
.size(); ++i
)
228 const int distCount
= sel_
[i
].posCount() / 2;
229 distances_
.setColumnCount(i
, distCount
);
230 xyz_
.setColumnCount(i
, distCount
* 3);
232 const double histogramMin
= (1.0 - lengthDev_
) * meanLength_
;
233 const double histogramMax
= (1.0 + lengthDev_
) * meanLength_
;
234 histogramModule_
->init(histogramFromRange(histogramMin
, histogramMax
)
235 .binWidth(binWidth_
).includeAll());
237 if (!fnAverage_
.empty())
239 AnalysisDataPlotModulePointer
plotm(
240 new AnalysisDataPlotModule(settings
.plotSettings()));
241 plotm
->setFileName(fnAverage_
);
242 plotm
->setTitle("Average distance");
243 plotm
->setXAxisIsTime();
244 plotm
->setYLabel("Distance (nm)");
245 for (size_t g
= 0; g
< sel_
.size(); ++g
)
247 plotm
->appendLegend(sel_
[g
].name());
249 averageModule_
->addModule(plotm
);
254 AnalysisDataPlotModulePointer
plotm(
255 new AnalysisDataPlotModule(settings
.plotSettings()));
256 plotm
->setFileName(fnAll_
);
257 plotm
->setTitle("Distance");
258 plotm
->setXAxisIsTime();
259 plotm
->setYLabel("Distance (nm)");
260 // TODO: Add legends? (there can be a massive amount of columns)
261 distances_
.addModule(plotm
);
266 AnalysisDataPlotModulePointer
plotm(
267 new AnalysisDataPlotModule(settings
.plotSettings()));
268 plotm
->setFileName(fnAll_
);
269 plotm
->setTitle("Distance");
270 plotm
->setXAxisIsTime();
271 plotm
->setYLabel("Distance (nm)");
272 // TODO: Add legends? (there can be a massive amount of columns)
273 xyz_
.addModule(plotm
);
276 if (!fnHistogram_
.empty())
278 AnalysisDataPlotModulePointer
plotm(
279 new AnalysisDataPlotModule(settings
.plotSettings()));
280 plotm
->setFileName(fnHistogram_
);
281 plotm
->setTitle("Distance histogram");
282 plotm
->setXLabel("Distance (nm)");
283 plotm
->setYLabel("Probability");
284 for (size_t g
= 0; g
< sel_
.size(); ++g
)
286 plotm
->appendLegend(sel_
[g
].name());
288 histogramModule_
->averager().addModule(plotm
);
291 if (!fnAllStats_
.empty())
293 AnalysisDataPlotModulePointer
plotm(
294 new AnalysisDataPlotModule(settings
.plotSettings()));
295 plotm
->setFileName(fnAllStats_
);
296 plotm
->setErrorsAsSeparateColumn(true);
297 plotm
->setTitle("Statistics for individual distances");
298 plotm
->setXLabel("Distance index");
299 plotm
->setYLabel("Average/standard deviation (nm)");
300 for (size_t g
= 0; g
< sel_
.size(); ++g
)
302 plotm
->appendLegend(std::string(sel_
[g
].name()) + " avg");
303 plotm
->appendLegend(std::string(sel_
[g
].name()) + " std.dev.");
305 // TODO: Consider whether this output format is the best possible.
306 allStatsModule_
->addModule(plotm
);
312 Distance::analyzeFrame(int frnr
, const t_trxframe
&fr
, t_pbc
*pbc
,
313 TrajectoryAnalysisModuleData
*pdata
)
315 AnalysisDataHandle distHandle
= pdata
->dataHandle(distances_
);
316 AnalysisDataHandle xyzHandle
= pdata
->dataHandle(xyz_
);
317 const SelectionList
&sel
= pdata
->parallelSelections(sel_
);
319 checkSelections(sel
);
321 distHandle
.startFrame(frnr
, fr
.time
);
322 xyzHandle
.startFrame(frnr
, fr
.time
);
323 for (size_t g
= 0; g
< sel
.size(); ++g
)
325 distHandle
.selectDataSet(g
);
326 xyzHandle
.selectDataSet(g
);
327 for (int i
= 0, n
= 0; i
< sel
[g
].posCount(); i
+= 2, ++n
)
329 const SelectionPosition
&p1
= sel
[g
].position(i
);
330 const SelectionPosition
&p2
= sel
[g
].position(i
+1);
334 pbc_dx(pbc
, p2
.x(), p1
.x(), dx
);
338 rvec_sub(p2
.x(), p1
.x(), dx
);
340 real dist
= norm(dx
);
341 bool bPresent
= p1
.selected() && p2
.selected();
342 distHandle
.setPoint(n
, dist
, bPresent
);
343 xyzHandle
.setPoints(n
*3, 3, dx
);
346 distHandle
.finishFrame();
347 xyzHandle
.finishFrame();
352 Distance::finishAnalysis(int /*nframes*/)
354 AbstractAverageHistogram
&averageHistogram
= histogramModule_
->averager();
355 averageHistogram
.normalizeProbability();
356 averageHistogram
.done();
361 Distance::writeOutput()
363 SelectionList::const_iterator sel
;
365 for (sel
= sel_
.begin(), index
= 0; sel
!= sel_
.end(); ++sel
, ++index
)
367 printf("%s:\n", sel
->name());
368 printf(" Number of samples: %d\n",
369 summaryStatsModule_
->sampleCount(index
, 0));
370 printf(" Average distance: %-6.3f nm\n",
371 summaryStatsModule_
->average(index
, 0));
372 printf(" Standard deviation: %-6.3f nm\n",
373 summaryStatsModule_
->standardDeviation(index
, 0));
379 const char DistanceInfo::name
[] = "distance";
380 const char DistanceInfo::shortDescription
[] =
381 "Calculate distances between pairs of positions";
383 TrajectoryAnalysisModulePointer
DistanceInfo::create()
385 return TrajectoryAnalysisModulePointer(new Distance
);
388 } // namespace analysismodules