2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 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::Angle.
39 * \author Teemu Murtola <teemu.murtola@gmail.com>
40 * \ingroup module_trajectoryanalysis
48 #include "gromacs/analysisdata/analysisdata.h"
49 #include "gromacs/analysisdata/modules/average.h"
50 #include "gromacs/analysisdata/modules/histogram.h"
51 #include "gromacs/analysisdata/modules/plot.h"
52 #include "gromacs/fileio/trx.h"
53 #include "gromacs/math/vec.h"
54 #include "gromacs/options/basicoptions.h"
55 #include "gromacs/options/filenameoption.h"
56 #include "gromacs/options/options.h"
57 #include "gromacs/pbcutil/pbc.h"
58 #include "gromacs/selection/selection.h"
59 #include "gromacs/selection/selectionoption.h"
60 #include "gromacs/trajectoryanalysis/analysissettings.h"
61 #include "gromacs/utility/arrayref.h"
62 #include "gromacs/utility/exceptions.h"
63 #include "gromacs/utility/gmxassert.h"
64 #include "gromacs/utility/stringutil.h"
69 namespace analysismodules
75 /********************************************************************
80 * Helper to encapsulate logic for looping over input selections.
82 * This class provides two-dimensional iteration:
83 * - Over _angle groups_, corresponding to an input selection. If the input
84 * selection list contains a single selection, that selection gets used
85 * for all angle groups.
86 * - Within a group, over _values_, each consisting of a fixed number of
87 * selection positions. If there is only a single value within a selection,
88 * that value is returned over and over again.
89 * This transparently provides the semantics of using a single selection/vector
90 * to compute angles against multiple selections/vectors as described in the
93 * This class isn't perferctly self-contained and requires the caller to know
94 * some of the internals to use it properly, but it serves its purpose for this
95 * single analysis tool by simplifying the loops.
96 * Some methods have also been tailored to allow the caller to use it a bit
99 * \ingroup module_trajectoryanalysis
101 class AnglePositionIterator
105 * Creates an iterator to loop over input selection positions.
107 * \param[in] selections List of selections.
108 * \param[in] posCountPerValue Number of selection positions that
109 * constitute a single value for the iteration.
111 * If \p selections is empty, and/or \p posCountPerValue is zero, the
112 * iterator can still be advanced and hasValue()/hasSingleValue()
113 * called, but values cannot be accessed.
115 AnglePositionIterator(const SelectionList &selections,
116 int posCountPerValue)
117 : selections_(selections), posCountPerValue_(posCountPerValue),
118 currentSelection_(0), nextPosition_(0)
122 //! Advances the iterator to the next group of angles.
125 if (selections_.size() > 1)
131 //! Advances the iterator to the next angle in the current group.
134 if (!hasSingleValue())
136 nextPosition_ += posCountPerValue_;
141 * Returns whether this iterator represents any values.
143 * If the return value is `false`, only nextGroup(), nextValue() and
144 * hasSingleValue() are allowed to be called.
146 bool hasValue() const
148 return !selections_.empty();
151 * Returns whether the current selection only contains a single value.
153 * Returns `false` if hasValue() returns false, which allows cutting
154 * some corners in consistency checks.
156 bool hasSingleValue() const
158 return hasValue() && currentSelection().posCount() == posCountPerValue_;
160 //! Returns whether the current selection is dynamic.
161 bool isDynamic() const
163 return currentSelection().isDynamic();
166 * Returns whether positions in the current value are either all
167 * selected or all unselected.
169 bool allValuesConsistentlySelected() const
171 if (posCountPerValue_ <= 1)
175 const bool bSelected = currentPosition(0).selected();
176 for (int i = 1; i < posCountPerValue_; ++i)
178 if (currentPosition(i).selected() != bSelected)
186 * Returns whether positions in the current value are selected.
188 * Only works reliably if allValuesConsistentlySelected() returns
191 bool currentValuesSelected() const
193 return selections_.empty() || currentPosition(0).selected();
196 //! Returns the currently active selection.
197 const Selection ¤tSelection() const
199 GMX_ASSERT(currentSelection_ < static_cast<int>(selections_.size()),
200 "Accessing an invalid selection");
201 return selections_[currentSelection_];
203 //! Returns the `i`th position for the current value.
204 SelectionPosition currentPosition(int i) const
206 return currentSelection().position(nextPosition_ + i);
209 * Extracts all coordinates corresponding to the current value.
211 * \param[out] x Array to which the positions are extracted.
213 * \p x should contain at minimum the number of positions per value
214 * passed to the constructor.
216 void getCurrentPositions(rvec x[]) const
218 GMX_ASSERT(posCountPerValue_ > 0,
219 "Accessing positions for an invalid angle type");
220 GMX_ASSERT(nextPosition_ + posCountPerValue_ <= currentSelection().posCount(),
221 "Accessing an invalid position");
222 for (int i = 0; i < posCountPerValue_; ++i)
224 copy_rvec(currentPosition(i).x(), x[i]);
229 const SelectionList &selections_;
230 const int posCountPerValue_;
231 int currentSelection_;
234 GMX_DISALLOW_COPY_AND_ASSIGN(AnglePositionIterator);
237 /********************************************************************
238 * Actual analysis module
241 class Angle : public TrajectoryAnalysisModule
247 virtual void initOptions(Options *options,
248 TrajectoryAnalysisSettings *settings);
249 virtual void optionsFinished(Options *options,
250 TrajectoryAnalysisSettings *settings);
251 virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
252 const TopologyInformation &top);
254 virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
255 TrajectoryAnalysisModuleData *pdata);
257 virtual void finishAnalysis(int nframes);
258 virtual void writeOutput();
261 void initFromSelections(const SelectionList &sel1,
262 const SelectionList &sel2);
263 void checkSelections(const SelectionList &sel1,
264 const SelectionList &sel2) const;
268 SelectionOptionInfo *sel1info_;
269 SelectionOptionInfo *sel2info_;
270 std::string fnAverage_;
272 std::string fnHistogram_;
278 AnalysisData angles_;
279 AnalysisDataFrameAverageModulePointer averageModule_;
280 AnalysisDataSimpleHistogramModulePointer histogramModule_;
282 std::vector<int> angleCount_;
285 // TODO: It is not possible to put rvec into a container.
286 std::vector<rvec *> vt0_;
288 // Copy and assign disallowed by base.
292 : TrajectoryAnalysisModule(AngleInfo::name, AngleInfo::shortDescription),
293 sel1info_(NULL), sel2info_(NULL), binWidth_(1.0), natoms1_(0), natoms2_(0)
295 averageModule_.reset(new AnalysisDataFrameAverageModule());
296 angles_.addModule(averageModule_);
297 histogramModule_.reset(new AnalysisDataSimpleHistogramModule());
298 angles_.addModule(histogramModule_);
300 registerAnalysisDataset(&angles_, "angle");
301 registerBasicDataset(averageModule_.get(), "average");
302 registerBasicDataset(&histogramModule_->averager(), "histogram");
308 for (size_t g = 0; g < vt0_.size(); ++g)
316 Angle::initOptions(Options *options, TrajectoryAnalysisSettings * /*settings*/)
318 static const char *const desc[] = {
319 "[THISMODULE] computes different types of angles between vectors.",
320 "It supports both vectors defined by two positions and normals of",
321 "planes defined by three positions.",
322 "The z axis or the local normal of a sphere can also be used as",
323 "one of the vectors.",
324 "There are also convenience options 'angle' and 'dihedral' for",
325 "calculating bond angles and dihedrals defined by three/four",
327 "The type of the angle is specified with [TT]-g1[tt] and [TT]-g2[tt].",
328 "If [TT]-g1[tt] is [TT]angle[tt] or [TT]dihedral[tt], [TT]-g2[tt]",
329 "should not be specified.",
330 "In this case, [TT]-group1[tt] should specify one or more selections,",
331 "and each should contain triplets or quartets of positions that define",
332 "the angles to be calculated.[PAR]",
333 "If [TT]-g1[tt] is [TT]vector[tt] or [TT]plane[tt], [TT]-group1[tt]",
334 "should specify selections that contain either pairs ([TT]vector[tt])",
335 "or triplets ([TT]plane[tt]) of positions. For vectors, the positions",
336 "set the endpoints of the vector, and for planes, the three positions",
337 "are used to calculate the normal of the plane. In both cases,",
338 "[TT]-g2[tt] specifies the other vector to use (see below).[PAR]",
339 "With [TT]-g2 vector[tt] or [TT]-g2 plane[tt], [TT]-group2[tt] should",
340 "specify another set of vectors. [TT]-group1[tt] and [TT]-group2[tt]",
341 "should specify the same number of selections. It is also allowed to",
342 "only have a single selection for one of the options, in which case",
343 "the same selection is used with each selection in the other group.",
344 "Similarly, for each selection in [TT]-group1[tt], the corresponding",
345 "selection in [TT]-group2[tt] should specify the same number of",
346 "vectors or a single vector. In the latter case, the angle is",
347 "calculated between that single vector and each vector from the other",
349 "With [TT]-g2 sphnorm[tt], each selection in [TT]-group2[tt] should",
350 "specify a single position that is the center of the sphere.",
351 "The second vector is calculated as the vector from the center to the",
352 "midpoint of the positions specified by [TT]-group1[tt].[PAR]",
353 "With [TT]-g2 z[tt], [TT]-group2[tt] is not necessary, and angles",
354 "between the first vectors and the positive Z axis are calculated.[PAR]",
355 "With [TT]-g2 t0[tt], [TT]-group2[tt] is not necessary, and angles",
356 "are calculated from the vectors as they are in the first frame.[PAR]",
357 "There are three options for output:",
358 "[TT]-oav[tt] writes an xvg file with the time and the average angle",
360 "[TT]-oall[tt] writes all the individual angles.",
361 "[TT]-oh[tt] writes a histogram of the angles. The bin width can be",
362 "set with [TT]-binw[tt].",
363 "For [TT]-oav[tt] and [TT]-oh[tt], separate average/histogram is",
364 "computed for each selection in [TT]-group1[tt]."
366 static const char *const cGroup1TypeEnum[] =
367 { "angle", "dihedral", "vector", "plane" };
368 static const char *const cGroup2TypeEnum[] =
369 { "none", "vector", "plane", "t0", "z", "sphnorm" };
371 options->setDescription(desc);
373 options->addOption(FileNameOption("oav").filetype(eftPlot).outputFile()
374 .store(&fnAverage_).defaultBasename("angaver")
375 .description("Average angles as a function of time"));
376 options->addOption(FileNameOption("oall").filetype(eftPlot).outputFile()
377 .store(&fnAll_).defaultBasename("angles")
378 .description("All angles as a function of time"));
379 options->addOption(FileNameOption("oh").filetype(eftPlot).outputFile()
380 .store(&fnHistogram_).defaultBasename("anghist")
381 .description("Histogram of the angles"));
383 options->addOption(StringOption("g1").enumValue(cGroup1TypeEnum)
384 .defaultEnumIndex(0).store(&g1type_)
385 .description("Type of analysis/first vector group"));
386 options->addOption(StringOption("g2").enumValue(cGroup2TypeEnum)
387 .defaultEnumIndex(0).store(&g2type_)
388 .description("Type of second vector group"));
389 options->addOption(DoubleOption("binw").store(&binWidth_)
390 .description("Binwidth for -oh in degrees"));
392 sel1info_ = options->addOption(SelectionOption("group1")
393 .required().dynamicMask().storeVector(&sel1_)
395 .description("First analysis/vector selection"));
396 sel2info_ = options->addOption(SelectionOption("group2")
397 .dynamicMask().storeVector(&sel2_)
399 .description("Second analysis/vector selection"));
404 Angle::optionsFinished(Options *options, TrajectoryAnalysisSettings * /* settings */)
406 const bool bSingle = (g1type_[0] == 'a' || g1type_[0] == 'd');
408 if (bSingle && g2type_[0] != 'n')
410 GMX_THROW(InconsistentInputError("Cannot use a second group (-g2) with "
411 "-g1 angle or dihedral"));
413 if (bSingle && options->isSet("group2"))
415 GMX_THROW(InconsistentInputError("Cannot provide a second selection "
416 "(-group2) with -g1 angle or dihedral"));
418 if (!bSingle && g2type_[0] == 'n')
420 GMX_THROW(InconsistentInputError("Should specify a second group (-g2) "
421 "if the first group is not an angle or a dihedral"));
424 // Set up the number of positions per angle.
427 case 'a': natoms1_ = 3; break;
428 case 'd': natoms1_ = 4; break;
429 case 'v': natoms1_ = 2; break;
430 case 'p': natoms1_ = 3; break;
432 GMX_THROW(InternalError("invalid -g1 value"));
436 case 'n': natoms2_ = 0; break;
437 case 'v': natoms2_ = 2; break;
438 case 'p': natoms2_ = 3; break;
439 case 't': natoms2_ = 0; break;
440 case 'z': natoms2_ = 0; break;
441 case 's': natoms2_ = 1; break;
443 GMX_THROW(InternalError("invalid -g2 value"));
445 if (natoms2_ == 0 && options->isSet("group2"))
447 GMX_THROW(InconsistentInputError("Cannot provide a second selection (-group2) with -g2 t0 or z"));
449 // TODO: If bSingle is not set, the second selection option should be
455 Angle::initFromSelections(const SelectionList &sel1,
456 const SelectionList &sel2)
458 const int angleGroups = std::max(sel1.size(), sel2.size());
459 const bool bHasSecondSelection = natoms2_ > 0;
461 if (bHasSecondSelection && sel1.size() != sel2.size()
462 && std::min(sel1.size(), sel2.size()) != 1)
464 GMX_THROW(InconsistentInputError(
465 "-group1 and -group2 should specify the same number of selections"));
468 AnglePositionIterator iter1(sel1, natoms1_);
469 AnglePositionIterator iter2(sel2, natoms2_);
470 for (int g = 0; g < angleGroups; ++g, iter1.nextGroup(), iter2.nextGroup())
472 const int posCount1 = iter1.currentSelection().posCount();
473 if (natoms1_ > 1 && posCount1 % natoms1_ != 0)
475 GMX_THROW(InconsistentInputError(formatString(
476 "Number of positions in selection %d in the first group not divisible by %d",
477 static_cast<int>(g + 1), natoms1_)));
479 const int angleCount1 = posCount1 / natoms1_;
480 int angleCount = angleCount1;
482 if (bHasSecondSelection)
484 const int posCount2 = iter2.currentSelection().posCount();
485 if (natoms2_ > 1 && posCount2 % natoms2_ != 0)
487 GMX_THROW(InconsistentInputError(formatString(
488 "Number of positions in selection %d in the second group not divisible by %d",
489 static_cast<int>(g + 1), natoms2_)));
491 if (g2type_[0] == 's' && posCount2 != 1)
493 GMX_THROW(InconsistentInputError(
494 "The second group should contain a single position with -g2 sphnorm"));
497 const int angleCount2 = posCount2 / natoms2_;
498 angleCount = std::max(angleCount1, angleCount2);
499 if (angleCount1 != angleCount2
500 && std::min(angleCount1, angleCount2) != 1)
502 GMX_THROW(InconsistentInputError(
503 "Number of vectors defined by the two groups are not the same"));
506 angleCount_.push_back(angleCount);
512 Angle::checkSelections(const SelectionList &sel1,
513 const SelectionList &sel2) const
515 AnglePositionIterator iter1(sel1, natoms1_);
516 AnglePositionIterator iter2(sel2, natoms2_);
517 for (size_t g = 0; g < angleCount_.size(); ++g, iter1.nextGroup(), iter2.nextGroup())
519 if (iter1.isDynamic() || (iter2.hasValue() && iter2.isDynamic()))
521 for (int n = 0; n < angleCount_[g]; ++n, iter1.nextValue(), iter2.nextValue())
524 if (!iter1.allValuesConsistentlySelected())
528 if (!iter2.allValuesConsistentlySelected())
532 if (angleCount_[g] > 1)
534 if (iter1.hasSingleValue() && !iter1.currentValuesSelected())
538 if (iter2.hasSingleValue() && !iter2.currentValuesSelected())
544 && (angleCount_[g] == 1
545 || (!iter1.hasSingleValue() && !iter2.hasSingleValue()))
546 && iter1.currentValuesSelected() != iter2.currentValuesSelected())
552 std::string message =
553 formatString("Dynamic selection %d does not select "
554 "a consistent set of angles over the frames",
555 static_cast<int>(g + 1));
556 GMX_THROW(InconsistentInputError(message));
565 Angle::initAnalysis(const TrajectoryAnalysisSettings &settings,
566 const TopologyInformation & /* top */)
568 initFromSelections(sel1_, sel2_);
570 // checkSelections() ensures that both selection lists have the same size.
571 angles_.setDataSetCount(angleCount_.size());
572 for (size_t i = 0; i < angleCount_.size(); ++i)
574 angles_.setColumnCount(i, angleCount_[i]);
576 double histogramMin = (g1type_ == "dihedral" ? -180.0 : 0);
577 histogramModule_->init(histogramFromRange(histogramMin, 180.0)
578 .binWidth(binWidth_).includeAll());
582 vt0_.resize(sel1_.size());
583 for (size_t g = 0; g < sel1_.size(); ++g)
585 vt0_[g] = new rvec[sel1_[g].posCount() / natoms1_];
589 if (!fnAverage_.empty())
591 AnalysisDataPlotModulePointer plotm(
592 new AnalysisDataPlotModule(settings.plotSettings()));
593 plotm->setFileName(fnAverage_);
594 plotm->setTitle("Average angle");
595 plotm->setXAxisIsTime();
596 plotm->setYLabel("Angle (degrees)");
597 // TODO: Consider adding information about the second selection,
598 // and/or a subtitle describing what kind of angle this is.
599 for (size_t g = 0; g < sel1_.size(); ++g)
601 plotm->appendLegend(sel1_[g].name());
603 averageModule_->addModule(plotm);
608 AnalysisDataPlotModulePointer plotm(
609 new AnalysisDataPlotModule(settings.plotSettings()));
610 plotm->setFileName(fnAll_);
611 plotm->setTitle("Angle");
612 plotm->setXAxisIsTime();
613 plotm->setYLabel("Angle (degrees)");
614 // TODO: Add legends? (there can be a massive amount of columns)
615 angles_.addModule(plotm);
618 if (!fnHistogram_.empty())
620 AnalysisDataPlotModulePointer plotm(
621 new AnalysisDataPlotModule(settings.plotSettings()));
622 plotm->setFileName(fnHistogram_);
623 plotm->setTitle("Angle histogram");
624 plotm->setXLabel("Angle (degrees)");
625 plotm->setYLabel("Probability");
626 // TODO: Consider adding information about the second selection,
627 // and/or a subtitle describing what kind of angle this is.
628 for (size_t g = 0; g < sel1_.size(); ++g)
630 plotm->appendLegend(sel1_[g].name());
632 histogramModule_->averager().addModule(plotm);
637 //! Helper method to calculate a vector from two or three positions.
639 calc_vec(int natoms, rvec x[], t_pbc *pbc, rvec xout, rvec cout)
646 pbc_dx(pbc, x[1], x[0], xout);
650 rvec_sub(x[1], x[0], xout);
652 svmul(0.5, xout, cout);
653 rvec_add(x[0], cout, cout);
660 pbc_dx(pbc, x[1], x[0], v1);
661 pbc_dx(pbc, x[2], x[0], v2);
665 rvec_sub(x[1], x[0], v1);
666 rvec_sub(x[2], x[0], v2);
669 rvec_add(x[0], x[1], cout);
670 rvec_add(cout, x[2], cout);
671 svmul(1.0/3.0, cout, cout);
675 GMX_RELEASE_ASSERT(false, "Incorrectly initialized number of atoms");
681 Angle::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
682 TrajectoryAnalysisModuleData *pdata)
684 AnalysisDataHandle dh = pdata->dataHandle(angles_);
685 const SelectionList &sel1 = pdata->parallelSelections(sel1_);
686 const SelectionList &sel2 = pdata->parallelSelections(sel2_);
688 checkSelections(sel1, sel2);
690 dh.startFrame(frnr, fr.time);
692 AnglePositionIterator iter1(sel1, natoms1_);
693 AnglePositionIterator iter2(sel2, natoms2_);
694 for (size_t g = 0; g < angleCount_.size(); ++g, iter1.nextGroup(), iter2.nextGroup())
706 copy_rvec(sel2_[g].position(0).x(), c2);
710 for (int n = 0; n < angleCount_[g]; ++n, iter1.nextValue(), iter2.nextValue())
714 // checkSelections() ensures that this reflects all the involved
716 const bool bPresent =
717 iter1.currentValuesSelected() && iter2.currentValuesSelected();
718 iter1.getCurrentPositions(x);
724 pbc_dx(pbc, x[0], x[1], v1);
725 pbc_dx(pbc, x[2], x[1], v2);
729 rvec_sub(x[0], x[1], v1);
730 rvec_sub(x[2], x[1], v2);
732 angle = gmx_angle(v1, v2);
739 pbc_dx(pbc, x[0], x[1], dx[0]);
740 pbc_dx(pbc, x[2], x[1], dx[1]);
741 pbc_dx(pbc, x[2], x[3], dx[2]);
745 rvec_sub(x[0], x[1], dx[0]);
746 rvec_sub(x[2], x[1], dx[1]);
747 rvec_sub(x[2], x[3], dx[2]);
749 cprod(dx[0], dx[1], v1);
750 cprod(dx[1], dx[2], v2);
751 angle = gmx_angle(v1, v2);
752 real ipr = iprod(dx[0], v2);
761 calc_vec(natoms1_, x, pbc, v1, c1);
766 iter2.getCurrentPositions(x);
767 calc_vec(natoms2_, x, pbc, v2, c2);
770 // FIXME: This is not parallelizable.
773 copy_rvec(v1, vt0_[g][n]);
775 copy_rvec(vt0_[g][n], v2);
778 c1[XX] = c1[YY] = 0.0;
783 pbc_dx(pbc, c1, c2, v2);
787 rvec_sub(c1, c2, v2);
791 GMX_THROW(InternalError("invalid -g2 value"));
793 angle = gmx_angle(v1, v2);
796 GMX_THROW(InternalError("invalid -g1 value"));
798 dh.setPoint(n, angle * RAD2DEG, bPresent);
806 Angle::finishAnalysis(int /*nframes*/)
808 AbstractAverageHistogram &averageHistogram = histogramModule_->averager();
809 averageHistogram.normalizeProbability();
810 averageHistogram.done();
821 const char AngleInfo::name[] = "gangle";
822 const char AngleInfo::shortDescription[] =
825 TrajectoryAnalysisModulePointer AngleInfo::create()
827 return TrajectoryAnalysisModulePointer(new Angle);
830 } // namespace analysismodules