2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2011,2012,2013, by the GROMACS development team, led by
5 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
6 * others, as listed in the AUTHORS file in the top-level source
7 * 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::Angle.
39 * \author Teemu Murtola <teemu.murtola@gmail.com>
40 * \ingroup module_trajectoryanalysis
47 #include "gromacs/legacyheaders/pbc.h"
48 #include "gromacs/legacyheaders/vec.h"
50 #include "gromacs/analysisdata/analysisdata.h"
51 #include "gromacs/analysisdata/modules/average.h"
52 #include "gromacs/analysisdata/modules/histogram.h"
53 #include "gromacs/analysisdata/modules/plot.h"
54 #include "gromacs/options/basicoptions.h"
55 #include "gromacs/options/filenameoption.h"
56 #include "gromacs/options/options.h"
57 #include "gromacs/selection/selection.h"
58 #include "gromacs/selection/selectionoption.h"
59 #include "gromacs/trajectoryanalysis/analysissettings.h"
60 #include "gromacs/utility/exceptions.h"
61 #include "gromacs/utility/gmxassert.h"
62 #include "gromacs/utility/stringutil.h"
67 namespace analysismodules
73 class Angle : public TrajectoryAnalysisModule
79 virtual void initOptions(Options *options,
80 TrajectoryAnalysisSettings *settings);
81 virtual void optionsFinished(Options *options,
82 TrajectoryAnalysisSettings *settings);
83 virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
84 const TopologyInformation &top);
86 virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
87 TrajectoryAnalysisModuleData *pdata);
89 virtual void finishAnalysis(int nframes);
90 virtual void writeOutput();
93 void checkSelections(const SelectionList &sel1,
94 const SelectionList &sel2) const;
98 SelectionOptionInfo *sel1info_;
99 SelectionOptionInfo *sel2info_;
100 std::string fnAverage_;
102 std::string fnHistogram_;
108 AnalysisData angles_;
109 AnalysisDataFrameAverageModulePointer averageModule_;
110 AnalysisDataSimpleHistogramModulePointer histogramModule_;
113 // TODO: It is not possible to put rvec into a container.
114 std::vector<rvec *> vt0_;
116 // Copy and assign disallowed by base.
120 : TrajectoryAnalysisModule(AngleInfo::name, AngleInfo::shortDescription),
121 sel1info_(NULL), sel2info_(NULL), binWidth_(1.0), natoms1_(0), natoms2_(0)
123 averageModule_.reset(new AnalysisDataFrameAverageModule());
124 angles_.addModule(averageModule_);
125 histogramModule_.reset(new AnalysisDataSimpleHistogramModule());
126 angles_.addModule(histogramModule_);
128 registerAnalysisDataset(&angles_, "angle");
129 registerBasicDataset(averageModule_.get(), "average");
130 registerBasicDataset(&histogramModule_->averager(), "histogram");
136 for (size_t g = 0; g < vt0_.size(); ++g)
144 Angle::initOptions(Options *options, TrajectoryAnalysisSettings * /*settings*/)
146 static const char *const desc[] = {
147 "g_angle computes different types of angles between vectors.",
148 "It supports both vectors defined by two positions and normals of",
149 "planes defined by three positions.",
150 "The z axis or the local normal of a sphere can also be used as",
151 "one of the vectors.",
152 "There are also convenience options 'angle' and 'dihedral' for",
153 "calculating bond angles and dihedrals defined by three/four",
155 "The type of the angle is specified with [TT]-g1[tt] and [TT]-g2[tt].",
156 "If [TT]-g1[tt] is [TT]angle[tt] or [TT]dihedral[tt], [TT]-g2[tt]",
157 "should not be specified.",
158 "In this case, [TT]-group1[tt] should specify one or more selections,",
159 "and each should contain triplets or quartets of positions that define",
160 "the angles to be calculated.[PAR]",
161 "If [TT]-g1[tt] is [TT]vector[tt] or [TT]plane[tt], [TT]-group1[tt]",
162 "should specify selections that contain either pairs ([TT]vector[tt])",
163 "or triplets ([TT]plane[tt]) of positions. For vectors, the positions",
164 "set the endpoints of the vector, and for planes, the three positions",
165 "are used to calculate the normal of the plane. In both cases,",
166 "[TT]-g2[tt] specifies the other vector to use (see below).[PAR]",
167 "With [TT]-g2 vector[tt] or [TT]-g2 plane[tt], [TT]-group2[tt] should",
168 "specify another set of vectors. [TT]-group1[tt] and [TT]-group2[tt]",
169 "should specify the same number of selections, and for each selection",
170 "in [TT]-group1[tt], the corresponding selection in [TT]-group2[tt]",
171 "should specify the same number of vectors.[PAR]",
172 "With [TT]-g2 sphnorm[tt], each selection in [TT]-group2[tt] should",
173 "specify a single position that is the center of the sphere.",
174 "The second vector is calculated as the vector from the center to the",
175 "midpoint of the positions specified by [TT]-group1[tt].[PAR]",
176 "With [TT]-g2 z[tt], [TT]-group2[tt] is not necessary, and angles",
177 "between the first vectors and the positive Z axis are calculated.[PAR]",
178 "With [TT]-g2 t0[tt], [TT]-group2[tt] is not necessary, and angles",
179 "are calculated from the vectors as they are in the first frame.[PAR]",
180 "There are three options for output:",
181 "[TT]-oav[tt] writes an xvgr file with the time and the average angle",
183 "[TT]-oall[tt] writes all the individual angles.",
184 "[TT]-oh[tt] writes a histogram of the angles. The bin width can be",
185 "set with [TT]-binw[tt].",
186 "For [TT]-oav[tt] and [TT]-oh[tt], separate average/histogram is",
187 "computed for each selection in [TT]-group1[tt]."
189 static const char *const cGroup1TypeEnum[] =
190 { "angle", "dihedral", "vector", "plane" };
191 static const char *const cGroup2TypeEnum[] =
192 { "none", "vector", "plane", "t0", "z", "sphnorm" };
194 options->setDescription(concatenateStrings(desc));
196 options->addOption(FileNameOption("oav").filetype(eftPlot).outputFile()
197 .store(&fnAverage_).defaultBasename("angaver")
198 .description("Average angles as a function of time"));
199 options->addOption(FileNameOption("oall").filetype(eftPlot).outputFile()
200 .store(&fnAll_).defaultBasename("angles")
201 .description("All angles as a function of time"));
202 options->addOption(FileNameOption("oh").filetype(eftPlot).outputFile()
203 .store(&fnHistogram_).defaultBasename("anghist")
204 .description("Histogram of the angles"));
206 options->addOption(StringOption("g1").enumValue(cGroup1TypeEnum)
207 .defaultEnumIndex(0).store(&g1type_)
208 .description("Type of analysis/first vector group"));
209 options->addOption(StringOption("g2").enumValue(cGroup2TypeEnum)
210 .defaultEnumIndex(0).store(&g2type_)
211 .description("Type of second vector group"));
212 options->addOption(DoubleOption("binw").store(&binWidth_)
213 .description("Binwidth for -oh in degrees"));
215 sel1info_ = options->addOption(SelectionOption("group1")
216 .required().dynamicMask().storeVector(&sel1_)
218 .description("First analysis/vector selection"));
219 sel2info_ = options->addOption(SelectionOption("group2")
220 .dynamicMask().storeVector(&sel2_)
222 .description("Second analysis/vector selection"));
227 Angle::optionsFinished(Options *options, TrajectoryAnalysisSettings * /* settings */)
229 bool bSingle = (g1type_[0] == 'a' || g1type_[0] == 'd');
231 if (bSingle && g2type_[0] != 'n')
233 GMX_THROW(InconsistentInputError("Cannot use a second group (-g2) with "
234 "-g1 angle or dihedral"));
236 if (bSingle && options->isSet("group2"))
238 GMX_THROW(InconsistentInputError("Cannot provide a second selection "
239 "(-group2) with -g1 angle or dihedral"));
241 if (!bSingle && g2type_[0] == 'n')
243 GMX_THROW(InconsistentInputError("Should specify a second group (-g2) "
244 "if the first group is not an angle or a dihedral"));
247 // Set up the number of positions per angle.
250 case 'a': natoms1_ = 3; break;
251 case 'd': natoms1_ = 4; break;
252 case 'v': natoms1_ = 2; break;
253 case 'p': natoms1_ = 3; break;
255 GMX_THROW(InternalError("invalid -g1 value"));
259 case 'n': natoms2_ = 0; break;
260 case 'v': natoms2_ = 2; break;
261 case 'p': natoms2_ = 3; break;
262 case 't': natoms2_ = 0; break;
263 case 'z': natoms2_ = 0; break;
264 case 's': natoms2_ = 1; break;
266 GMX_THROW(InternalError("invalid -g2 value"));
268 if (natoms2_ == 0 && options->isSet("group2"))
270 GMX_THROW(InconsistentInputError("Cannot provide a second selection (-group2) with -g2 t0 or z"));
276 Angle::checkSelections(const SelectionList &sel1,
277 const SelectionList &sel2) const
279 if (natoms2_ > 0 && sel1.size() != sel2.size())
281 GMX_THROW(InconsistentInputError(
282 "-group1 and -group2 should specify the same number of selections"));
285 for (size_t g = 0; g < sel1.size(); ++g)
287 const int na1 = sel1[g].posCount();
288 const int na2 = (natoms2_ > 0) ? sel2[g].posCount() : 0;
289 if (natoms1_ > 1 && na1 % natoms1_ != 0)
291 GMX_THROW(InconsistentInputError(formatString(
292 "Number of positions in selection %d in the first group not divisible by %d",
293 static_cast<int>(g + 1), natoms1_)));
295 if (natoms2_ > 1 && na2 % natoms2_ != 0)
297 GMX_THROW(InconsistentInputError(formatString(
298 "Number of positions in selection %d in the second group not divisible by %d",
299 static_cast<int>(g + 1), natoms2_)));
301 if (natoms1_ > 0 && natoms2_ > 1 && na1 / natoms1_ != na2 / natoms2_)
303 GMX_THROW(InconsistentInputError(
304 "Number of vectors defined by the two groups are not the same"));
306 if (g2type_[0] == 's' && sel2[g].posCount() != 1)
308 GMX_THROW(InconsistentInputError(
309 "The second group should contain a single position with -g2 sphnorm"));
311 if (sel1[g].isDynamic() || (natoms2_ > 0 && sel2[g].isDynamic()))
313 for (int i = 0, j = 0; i < na1; i += natoms1_, j += natoms2_)
315 const bool bSelected = sel1[g].position(i).selected();
317 for (int k = 1; k < natoms1_ && bOk; ++k)
319 bOk = (sel1[g].position(i+k).selected() == bSelected);
321 for (int k = 1; k < natoms2_ && bOk; ++k)
323 bOk = (sel2[g].position(j+k).selected() == bSelected);
327 std::string message =
328 formatString("Dynamic selection %d does not select "
329 "a consistent set of angles over the frames",
330 static_cast<int>(g + 1));
331 GMX_THROW(InconsistentInputError(message));
340 Angle::initAnalysis(const TrajectoryAnalysisSettings &settings,
341 const TopologyInformation & /* top */)
343 checkSelections(sel1_, sel2_);
345 // checkSelections() ensures that both selection lists have the same size.
346 angles_.setDataSetCount(sel1_.size());
347 for (size_t i = 0; i < sel1_.size(); ++i)
349 angles_.setColumnCount(i, sel1_[i].posCount() / natoms1_);
351 double histogramMin = (g1type_ == "dihedral" ? -180.0 : 0);
352 histogramModule_->init(histogramFromRange(histogramMin, 180.0)
353 .binWidth(binWidth_).includeAll());
357 vt0_.resize(sel1_.size());
358 for (size_t g = 0; g < sel1_.size(); ++g)
360 vt0_[g] = new rvec[sel1_[g].posCount() / natoms1_];
364 if (!fnAverage_.empty())
366 AnalysisDataPlotModulePointer plotm(
367 new AnalysisDataPlotModule(settings.plotSettings()));
368 plotm->setFileName(fnAverage_);
369 plotm->setTitle("Average angle");
370 plotm->setXAxisIsTime();
371 plotm->setYLabel("Angle (degrees)");
372 // TODO: Consider adding information about the second selection,
373 // and/or a subtitle describing what kind of angle this is.
374 for (size_t g = 0; g < sel1_.size(); ++g)
376 plotm->appendLegend(sel1_[g].name());
378 averageModule_->addModule(plotm);
383 AnalysisDataPlotModulePointer plotm(
384 new AnalysisDataPlotModule(settings.plotSettings()));
385 plotm->setFileName(fnAll_);
386 plotm->setTitle("Angle");
387 plotm->setXAxisIsTime();
388 plotm->setYLabel("Angle (degrees)");
389 // TODO: Add legends? (there can be a massive amount of columns)
390 angles_.addModule(plotm);
393 if (!fnHistogram_.empty())
395 AnalysisDataPlotModulePointer plotm(
396 new AnalysisDataPlotModule(settings.plotSettings()));
397 plotm->setFileName(fnHistogram_);
398 plotm->setTitle("Angle histogram");
399 plotm->setXLabel("Angle (degrees)");
400 plotm->setYLabel("Probability");
401 // TODO: Consider adding information about the second selection,
402 // and/or a subtitle describing what kind of angle this is.
403 for (size_t g = 0; g < sel1_.size(); ++g)
405 plotm->appendLegend(sel1_[g].name());
407 histogramModule_->averager().addModule(plotm);
412 //! Helper method to process selections into an array of coordinates.
414 copy_pos(const SelectionList &sel, int natoms, int g, int first, rvec x[])
416 for (int k = 0; k < natoms; ++k)
418 copy_rvec(sel[g].position(first + k).x(), x[k]);
423 //! Helper method to calculate a vector from two or three positions.
425 calc_vec(int natoms, rvec x[], t_pbc *pbc, rvec xout, rvec cout)
432 pbc_dx(pbc, x[1], x[0], xout);
436 rvec_sub(x[1], x[0], xout);
438 svmul(0.5, xout, cout);
439 rvec_add(x[0], cout, cout);
446 pbc_dx(pbc, x[1], x[0], v1);
447 pbc_dx(pbc, x[2], x[0], v2);
451 rvec_sub(x[1], x[0], v1);
452 rvec_sub(x[2], x[0], v2);
455 rvec_add(x[0], x[1], cout);
456 rvec_add(cout, x[2], cout);
457 svmul(1.0/3.0, cout, cout);
461 GMX_RELEASE_ASSERT(false, "Incorrectly initialized number of atoms");
467 Angle::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
468 TrajectoryAnalysisModuleData *pdata)
470 AnalysisDataHandle dh = pdata->dataHandle(angles_);
471 const SelectionList &sel1 = pdata->parallelSelections(sel1_);
472 const SelectionList &sel2 = pdata->parallelSelections(sel2_);
474 checkSelections(sel1, sel2);
476 dh.startFrame(frnr, fr.time);
478 for (size_t g = 0; g < sel1_.size(); ++g)
490 copy_rvec(sel2_[g].position(0).x(), c2);
494 for (int i = 0, j = 0, n = 0;
495 i < sel1[g].posCount();
496 i += natoms1_, j += natoms2_, ++n)
500 // checkSelections() ensures that this reflects all the involved
502 bool bPresent = sel1[g].position(i).selected();
503 copy_pos(sel1, natoms1_, g, i, x);
509 pbc_dx(pbc, x[0], x[1], v1);
510 pbc_dx(pbc, x[2], x[1], v2);
514 rvec_sub(x[0], x[1], v1);
515 rvec_sub(x[2], x[1], v2);
517 angle = gmx_angle(v1, v2);
524 pbc_dx(pbc, x[0], x[1], dx[0]);
525 pbc_dx(pbc, x[2], x[1], dx[1]);
526 pbc_dx(pbc, x[2], x[3], dx[2]);
530 rvec_sub(x[0], x[1], dx[0]);
531 rvec_sub(x[2], x[1], dx[1]);
532 rvec_sub(x[2], x[3], dx[2]);
534 cprod(dx[0], dx[1], v1);
535 cprod(dx[1], dx[2], v2);
536 angle = gmx_angle(v1, v2);
537 real ipr = iprod(dx[0], v2);
546 calc_vec(natoms1_, x, pbc, v1, c1);
551 copy_pos(sel2, natoms2_, 0, j, x);
552 calc_vec(natoms2_, x, pbc, v2, c2);
555 // FIXME: This is not parallelizable.
558 copy_rvec(v1, vt0_[g][n]);
560 copy_rvec(vt0_[g][n], v2);
563 c1[XX] = c1[YY] = 0.0;
568 pbc_dx(pbc, c1, c2, v2);
572 rvec_sub(c1, c2, v2);
576 GMX_THROW(InternalError("invalid -g2 value"));
578 angle = gmx_angle(v1, v2);
581 GMX_THROW(InternalError("invalid -g1 value"));
583 dh.setPoint(n, angle * RAD2DEG, bPresent);
591 Angle::finishAnalysis(int /*nframes*/)
593 AbstractAverageHistogram &averageHistogram = histogramModule_->averager();
594 averageHistogram.normalizeProbability();
595 averageHistogram.done();
606 const char AngleInfo::name[] = "gangle";
607 const char AngleInfo::shortDescription[] =
610 TrajectoryAnalysisModulePointer AngleInfo::create()
612 return TrajectoryAnalysisModulePointer(new Angle);
615 } // namespace analysismodules