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
48 #include "gromacs/legacyheaders/pbc.h"
49 #include "gromacs/legacyheaders/vec.h"
51 #include "gromacs/analysisdata/analysisdata.h"
52 #include "gromacs/analysisdata/modules/average.h"
53 #include "gromacs/analysisdata/modules/histogram.h"
54 #include "gromacs/analysisdata/modules/plot.h"
55 #include "gromacs/options/basicoptions.h"
56 #include "gromacs/options/filenameoption.h"
57 #include "gromacs/options/options.h"
58 #include "gromacs/selection/selection.h"
59 #include "gromacs/selection/selectionoption.h"
60 #include "gromacs/trajectoryanalysis/analysissettings.h"
61 #include "gromacs/utility/exceptions.h"
62 #include "gromacs/utility/gmxassert.h"
63 #include "gromacs/utility/stringutil.h"
68 namespace analysismodules
74 /********************************************************************
79 * Helper to encapsulate logic for looping over input selections.
81 * This class provides two-dimensional iteration:
82 * - Over _angle groups_, corresponding to an input selection. If the input
83 * selection list contains a single selection, that selection gets used
84 * for all angle groups.
85 * - Within a group, over _values_, each consisting of a fixed number of
86 * selection positions. If there is only a single value within a selection,
87 * that value is returned over and over again.
88 * This transparently provides the semantics of using a single selection/vector
89 * to compute angles against multiple selections/vectors as described in the
92 * This class isn't perferctly self-contained and requires the caller to know
93 * some of the internals to use it properly, but it serves its purpose for this
94 * single analysis tool by simplifying the loops.
95 * Some methods have also been tailored to allow the caller to use it a bit
98 * \ingroup module_trajectoryanalysis
100 class AnglePositionIterator
104 * Creates an iterator to loop over input selection positions.
106 * \param[in] selections List of selections.
107 * \param[in] posCountPerValue Number of selection positions that
108 * constitute a single value for the iteration.
110 * If \p selections is empty, and/or \p posCountPerValue is zero, the
111 * iterator can still be advanced and hasValue()/hasSingleValue()
112 * called, but values cannot be accessed.
114 AnglePositionIterator(const SelectionList &selections,
115 int posCountPerValue)
116 : selections_(selections), posCountPerValue_(posCountPerValue),
117 currentSelection_(0), nextPosition_(0)
121 //! Advances the iterator to the next group of angles.
124 if (selections_.size() > 1)
130 //! Advances the iterator to the next angle in the current group.
133 if (!hasSingleValue())
135 nextPosition_ += posCountPerValue_;
140 * Returns whether this iterator represents any values.
142 * If the return value is `false`, only nextGroup(), nextValue() and
143 * hasSingleValue() are allowed to be called.
145 bool hasValue() const
147 return !selections_.empty();
150 * Returns whether the current selection only contains a single value.
152 * Returns `false` if hasValue() returns false, which allows cutting
153 * some corners in consistency checks.
155 bool hasSingleValue() const
157 return hasValue() && currentSelection().posCount() == posCountPerValue_;
159 //! Returns whether the current selection is dynamic.
160 bool isDynamic() const
162 return currentSelection().isDynamic();
165 * Returns whether positions in the current value are either all
166 * selected or all unselected.
168 bool allValuesConsistentlySelected() const
170 if (posCountPerValue_ <= 1)
174 const bool bSelected = currentPosition(0).selected();
175 for (int i = 1; i < posCountPerValue_; ++i)
177 if (currentPosition(i).selected() != bSelected)
185 * Returns whether positions in the current value are selected.
187 * Only works reliably if allValuesConsistentlySelected() returns
190 bool currentValuesSelected() const
192 return selections_.empty() || currentPosition(0).selected();
195 //! Returns the currently active selection.
196 const Selection ¤tSelection() const
198 GMX_ASSERT(currentSelection_ < static_cast<int>(selections_.size()),
199 "Accessing an invalid selection");
200 return selections_[currentSelection_];
202 //! Returns the `i`th position for the current value.
203 SelectionPosition currentPosition(int i) const
205 return currentSelection().position(nextPosition_ + i);
208 * Extracts all coordinates corresponding to the current value.
210 * \param[out] x Array to which the positions are extracted.
212 * \p x should contain at minimum the number of positions per value
213 * passed to the constructor.
215 void getCurrentPositions(rvec x[]) const
217 GMX_ASSERT(posCountPerValue_ > 0,
218 "Accessing positions for an invalid angle type");
219 GMX_ASSERT(nextPosition_ + posCountPerValue_ <= currentSelection().posCount(),
220 "Accessing an invalid position");
221 for (int i = 0; i < posCountPerValue_; ++i)
223 copy_rvec(currentPosition(i).x(), x[i]);
228 const SelectionList &selections_;
229 const int posCountPerValue_;
230 int currentSelection_;
233 GMX_DISALLOW_COPY_AND_ASSIGN(AnglePositionIterator);
236 /********************************************************************
237 * Actual analysis module
240 class Angle : public TrajectoryAnalysisModule
246 virtual void initOptions(Options *options,
247 TrajectoryAnalysisSettings *settings);
248 virtual void optionsFinished(Options *options,
249 TrajectoryAnalysisSettings *settings);
250 virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
251 const TopologyInformation &top);
253 virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
254 TrajectoryAnalysisModuleData *pdata);
256 virtual void finishAnalysis(int nframes);
257 virtual void writeOutput();
260 void initFromSelections(const SelectionList &sel1,
261 const SelectionList &sel2);
262 void checkSelections(const SelectionList &sel1,
263 const SelectionList &sel2) const;
267 SelectionOptionInfo *sel1info_;
268 SelectionOptionInfo *sel2info_;
269 std::string fnAverage_;
271 std::string fnHistogram_;
277 AnalysisData angles_;
278 AnalysisDataFrameAverageModulePointer averageModule_;
279 AnalysisDataSimpleHistogramModulePointer histogramModule_;
281 std::vector<int> angleCount_;
284 // TODO: It is not possible to put rvec into a container.
285 std::vector<rvec *> vt0_;
287 // Copy and assign disallowed by base.
291 : TrajectoryAnalysisModule(AngleInfo::name, AngleInfo::shortDescription),
292 sel1info_(NULL), sel2info_(NULL), binWidth_(1.0), natoms1_(0), natoms2_(0)
294 averageModule_.reset(new AnalysisDataFrameAverageModule());
295 angles_.addModule(averageModule_);
296 histogramModule_.reset(new AnalysisDataSimpleHistogramModule());
297 angles_.addModule(histogramModule_);
299 registerAnalysisDataset(&angles_, "angle");
300 registerBasicDataset(averageModule_.get(), "average");
301 registerBasicDataset(&histogramModule_->averager(), "histogram");
307 for (size_t g = 0; g < vt0_.size(); ++g)
315 Angle::initOptions(Options *options, TrajectoryAnalysisSettings * /*settings*/)
317 static const char *const desc[] = {
318 "[TT]gmx gangle[tt] computes different types of angles between vectors.",
319 "It supports both vectors defined by two positions and normals of",
320 "planes defined by three positions.",
321 "The z axis or the local normal of a sphere can also be used as",
322 "one of the vectors.",
323 "There are also convenience options 'angle' and 'dihedral' for",
324 "calculating bond angles and dihedrals defined by three/four",
326 "The type of the angle is specified with [TT]-g1[tt] and [TT]-g2[tt].",
327 "If [TT]-g1[tt] is [TT]angle[tt] or [TT]dihedral[tt], [TT]-g2[tt]",
328 "should not be specified.",
329 "In this case, [TT]-group1[tt] should specify one or more selections,",
330 "and each should contain triplets or quartets of positions that define",
331 "the angles to be calculated.[PAR]",
332 "If [TT]-g1[tt] is [TT]vector[tt] or [TT]plane[tt], [TT]-group1[tt]",
333 "should specify selections that contain either pairs ([TT]vector[tt])",
334 "or triplets ([TT]plane[tt]) of positions. For vectors, the positions",
335 "set the endpoints of the vector, and for planes, the three positions",
336 "are used to calculate the normal of the plane. In both cases,",
337 "[TT]-g2[tt] specifies the other vector to use (see below).[PAR]",
338 "With [TT]-g2 vector[tt] or [TT]-g2 plane[tt], [TT]-group2[tt] should",
339 "specify another set of vectors. [TT]-group1[tt] and [TT]-group2[tt]",
340 "should specify the same number of selections. It is also allowed to",
341 "only have a single selection for one of the options, in which case",
342 "the same selection is used with each selection in the other group.",
343 "Similarly, for each selection in [TT]-group1[tt], the corresponding",
344 "selection in [TT]-group2[tt] should specify the same number of",
345 "vectors or a single vector. In the latter case, the angle is",
346 "calculated between that single vector and each vector from the other",
348 "With [TT]-g2 sphnorm[tt], each selection in [TT]-group2[tt] should",
349 "specify a single position that is the center of the sphere.",
350 "The second vector is calculated as the vector from the center to the",
351 "midpoint of the positions specified by [TT]-group1[tt].[PAR]",
352 "With [TT]-g2 z[tt], [TT]-group2[tt] is not necessary, and angles",
353 "between the first vectors and the positive Z axis are calculated.[PAR]",
354 "With [TT]-g2 t0[tt], [TT]-group2[tt] is not necessary, and angles",
355 "are calculated from the vectors as they are in the first frame.[PAR]",
356 "There are three options for output:",
357 "[TT]-oav[tt] writes an xvgr file with the time and the average angle",
359 "[TT]-oall[tt] writes all the individual angles.",
360 "[TT]-oh[tt] writes a histogram of the angles. The bin width can be",
361 "set with [TT]-binw[tt].",
362 "For [TT]-oav[tt] and [TT]-oh[tt], separate average/histogram is",
363 "computed for each selection in [TT]-group1[tt]."
365 static const char *const cGroup1TypeEnum[] =
366 { "angle", "dihedral", "vector", "plane" };
367 static const char *const cGroup2TypeEnum[] =
368 { "none", "vector", "plane", "t0", "z", "sphnorm" };
370 options->setDescription(concatenateStrings(desc));
372 options->addOption(FileNameOption("oav").filetype(eftPlot).outputFile()
373 .store(&fnAverage_).defaultBasename("angaver")
374 .description("Average angles as a function of time"));
375 options->addOption(FileNameOption("oall").filetype(eftPlot).outputFile()
376 .store(&fnAll_).defaultBasename("angles")
377 .description("All angles as a function of time"));
378 options->addOption(FileNameOption("oh").filetype(eftPlot).outputFile()
379 .store(&fnHistogram_).defaultBasename("anghist")
380 .description("Histogram of the angles"));
382 options->addOption(StringOption("g1").enumValue(cGroup1TypeEnum)
383 .defaultEnumIndex(0).store(&g1type_)
384 .description("Type of analysis/first vector group"));
385 options->addOption(StringOption("g2").enumValue(cGroup2TypeEnum)
386 .defaultEnumIndex(0).store(&g2type_)
387 .description("Type of second vector group"));
388 options->addOption(DoubleOption("binw").store(&binWidth_)
389 .description("Binwidth for -oh in degrees"));
391 sel1info_ = options->addOption(SelectionOption("group1")
392 .required().dynamicMask().storeVector(&sel1_)
394 .description("First analysis/vector selection"));
395 sel2info_ = options->addOption(SelectionOption("group2")
396 .dynamicMask().storeVector(&sel2_)
398 .description("Second analysis/vector selection"));
403 Angle::optionsFinished(Options *options, TrajectoryAnalysisSettings * /* settings */)
405 const bool bSingle = (g1type_[0] == 'a' || g1type_[0] == 'd');
407 if (bSingle && g2type_[0] != 'n')
409 GMX_THROW(InconsistentInputError("Cannot use a second group (-g2) with "
410 "-g1 angle or dihedral"));
412 if (bSingle && options->isSet("group2"))
414 GMX_THROW(InconsistentInputError("Cannot provide a second selection "
415 "(-group2) with -g1 angle or dihedral"));
417 if (!bSingle && g2type_[0] == 'n')
419 GMX_THROW(InconsistentInputError("Should specify a second group (-g2) "
420 "if the first group is not an angle or a dihedral"));
423 // Set up the number of positions per angle.
426 case 'a': natoms1_ = 3; break;
427 case 'd': natoms1_ = 4; break;
428 case 'v': natoms1_ = 2; break;
429 case 'p': natoms1_ = 3; break;
431 GMX_THROW(InternalError("invalid -g1 value"));
435 case 'n': natoms2_ = 0; break;
436 case 'v': natoms2_ = 2; break;
437 case 'p': natoms2_ = 3; break;
438 case 't': natoms2_ = 0; break;
439 case 'z': natoms2_ = 0; break;
440 case 's': natoms2_ = 1; break;
442 GMX_THROW(InternalError("invalid -g2 value"));
444 if (natoms2_ == 0 && options->isSet("group2"))
446 GMX_THROW(InconsistentInputError("Cannot provide a second selection (-group2) with -g2 t0 or z"));
448 // TODO: If bSingle is not set, the second selection option should be
454 Angle::initFromSelections(const SelectionList &sel1,
455 const SelectionList &sel2)
457 const int angleGroups = std::max(sel1.size(), sel2.size());
458 const bool bHasSecondSelection = natoms2_ > 0;
460 if (bHasSecondSelection && sel1.size() != sel2.size()
461 && std::min(sel1.size(), sel2.size()) != 1)
463 GMX_THROW(InconsistentInputError(
464 "-group1 and -group2 should specify the same number of selections"));
467 AnglePositionIterator iter1(sel1, natoms1_);
468 AnglePositionIterator iter2(sel2, natoms2_);
469 for (int g = 0; g < angleGroups; ++g, iter1.nextGroup(), iter2.nextGroup())
471 const int posCount1 = iter1.currentSelection().posCount();
472 if (natoms1_ > 1 && posCount1 % natoms1_ != 0)
474 GMX_THROW(InconsistentInputError(formatString(
475 "Number of positions in selection %d in the first group not divisible by %d",
476 static_cast<int>(g + 1), natoms1_)));
478 const int angleCount1 = posCount1 / natoms1_;
479 int angleCount = angleCount1;
481 if (bHasSecondSelection)
483 const int posCount2 = iter2.currentSelection().posCount();
484 if (natoms2_ > 1 && posCount2 % natoms2_ != 0)
486 GMX_THROW(InconsistentInputError(formatString(
487 "Number of positions in selection %d in the second group not divisible by %d",
488 static_cast<int>(g + 1), natoms2_)));
490 if (g2type_[0] == 's' && posCount2 != 1)
492 GMX_THROW(InconsistentInputError(
493 "The second group should contain a single position with -g2 sphnorm"));
496 const int angleCount2 = posCount2 / natoms2_;
497 angleCount = std::max(angleCount1, angleCount2);
498 if (angleCount1 != angleCount2
499 && std::min(angleCount1, angleCount2) != 1)
501 GMX_THROW(InconsistentInputError(
502 "Number of vectors defined by the two groups are not the same"));
505 angleCount_.push_back(angleCount);
511 Angle::checkSelections(const SelectionList &sel1,
512 const SelectionList &sel2) const
514 AnglePositionIterator iter1(sel1, natoms1_);
515 AnglePositionIterator iter2(sel2, natoms2_);
516 for (size_t g = 0; g < angleCount_.size(); ++g, iter1.nextGroup(), iter2.nextGroup())
518 if (iter1.isDynamic() || (iter2.hasValue() && iter2.isDynamic()))
520 for (int n = 0; n < angleCount_[g]; ++n, iter1.nextValue(), iter2.nextValue())
523 if (!iter1.allValuesConsistentlySelected())
527 if (!iter2.allValuesConsistentlySelected())
531 if (angleCount_[g] > 1)
533 if (iter1.hasSingleValue() && !iter1.currentValuesSelected())
537 if (iter2.hasSingleValue() && !iter2.currentValuesSelected())
543 && (angleCount_[g] == 1
544 || (!iter1.hasSingleValue() && !iter2.hasSingleValue()))
545 && iter1.currentValuesSelected() != iter2.currentValuesSelected())
551 std::string message =
552 formatString("Dynamic selection %d does not select "
553 "a consistent set of angles over the frames",
554 static_cast<int>(g + 1));
555 GMX_THROW(InconsistentInputError(message));
564 Angle::initAnalysis(const TrajectoryAnalysisSettings &settings,
565 const TopologyInformation & /* top */)
567 initFromSelections(sel1_, sel2_);
569 // checkSelections() ensures that both selection lists have the same size.
570 angles_.setDataSetCount(angleCount_.size());
571 for (size_t i = 0; i < angleCount_.size(); ++i)
573 angles_.setColumnCount(i, angleCount_[i]);
575 double histogramMin = (g1type_ == "dihedral" ? -180.0 : 0);
576 histogramModule_->init(histogramFromRange(histogramMin, 180.0)
577 .binWidth(binWidth_).includeAll());
581 vt0_.resize(sel1_.size());
582 for (size_t g = 0; g < sel1_.size(); ++g)
584 vt0_[g] = new rvec[sel1_[g].posCount() / natoms1_];
588 if (!fnAverage_.empty())
590 AnalysisDataPlotModulePointer plotm(
591 new AnalysisDataPlotModule(settings.plotSettings()));
592 plotm->setFileName(fnAverage_);
593 plotm->setTitle("Average angle");
594 plotm->setXAxisIsTime();
595 plotm->setYLabel("Angle (degrees)");
596 // TODO: Consider adding information about the second selection,
597 // and/or a subtitle describing what kind of angle this is.
598 for (size_t g = 0; g < sel1_.size(); ++g)
600 plotm->appendLegend(sel1_[g].name());
602 averageModule_->addModule(plotm);
607 AnalysisDataPlotModulePointer plotm(
608 new AnalysisDataPlotModule(settings.plotSettings()));
609 plotm->setFileName(fnAll_);
610 plotm->setTitle("Angle");
611 plotm->setXAxisIsTime();
612 plotm->setYLabel("Angle (degrees)");
613 // TODO: Add legends? (there can be a massive amount of columns)
614 angles_.addModule(plotm);
617 if (!fnHistogram_.empty())
619 AnalysisDataPlotModulePointer plotm(
620 new AnalysisDataPlotModule(settings.plotSettings()));
621 plotm->setFileName(fnHistogram_);
622 plotm->setTitle("Angle histogram");
623 plotm->setXLabel("Angle (degrees)");
624 plotm->setYLabel("Probability");
625 // TODO: Consider adding information about the second selection,
626 // and/or a subtitle describing what kind of angle this is.
627 for (size_t g = 0; g < sel1_.size(); ++g)
629 plotm->appendLegend(sel1_[g].name());
631 histogramModule_->averager().addModule(plotm);
636 //! Helper method to calculate a vector from two or three positions.
638 calc_vec(int natoms, rvec x[], t_pbc *pbc, rvec xout, rvec cout)
645 pbc_dx(pbc, x[1], x[0], xout);
649 rvec_sub(x[1], x[0], xout);
651 svmul(0.5, xout, cout);
652 rvec_add(x[0], cout, cout);
659 pbc_dx(pbc, x[1], x[0], v1);
660 pbc_dx(pbc, x[2], x[0], v2);
664 rvec_sub(x[1], x[0], v1);
665 rvec_sub(x[2], x[0], v2);
668 rvec_add(x[0], x[1], cout);
669 rvec_add(cout, x[2], cout);
670 svmul(1.0/3.0, cout, cout);
674 GMX_RELEASE_ASSERT(false, "Incorrectly initialized number of atoms");
680 Angle::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
681 TrajectoryAnalysisModuleData *pdata)
683 AnalysisDataHandle dh = pdata->dataHandle(angles_);
684 const SelectionList &sel1 = pdata->parallelSelections(sel1_);
685 const SelectionList &sel2 = pdata->parallelSelections(sel2_);
687 checkSelections(sel1, sel2);
689 dh.startFrame(frnr, fr.time);
691 AnglePositionIterator iter1(sel1, natoms1_);
692 AnglePositionIterator iter2(sel2, natoms2_);
693 for (size_t g = 0; g < angleCount_.size(); ++g, iter1.nextGroup(), iter2.nextGroup())
705 copy_rvec(sel2_[g].position(0).x(), c2);
709 for (int n = 0; n < angleCount_[g]; ++n, iter1.nextValue(), iter2.nextValue())
713 // checkSelections() ensures that this reflects all the involved
715 const bool bPresent =
716 iter1.currentValuesSelected() && iter2.currentValuesSelected();
717 iter1.getCurrentPositions(x);
723 pbc_dx(pbc, x[0], x[1], v1);
724 pbc_dx(pbc, x[2], x[1], v2);
728 rvec_sub(x[0], x[1], v1);
729 rvec_sub(x[2], x[1], v2);
731 angle = gmx_angle(v1, v2);
738 pbc_dx(pbc, x[0], x[1], dx[0]);
739 pbc_dx(pbc, x[2], x[1], dx[1]);
740 pbc_dx(pbc, x[2], x[3], dx[2]);
744 rvec_sub(x[0], x[1], dx[0]);
745 rvec_sub(x[2], x[1], dx[1]);
746 rvec_sub(x[2], x[3], dx[2]);
748 cprod(dx[0], dx[1], v1);
749 cprod(dx[1], dx[2], v2);
750 angle = gmx_angle(v1, v2);
751 real ipr = iprod(dx[0], v2);
760 calc_vec(natoms1_, x, pbc, v1, c1);
765 iter2.getCurrentPositions(x);
766 calc_vec(natoms2_, x, pbc, v2, c2);
769 // FIXME: This is not parallelizable.
772 copy_rvec(v1, vt0_[g][n]);
774 copy_rvec(vt0_[g][n], v2);
777 c1[XX] = c1[YY] = 0.0;
782 pbc_dx(pbc, c1, c2, v2);
786 rvec_sub(c1, c2, v2);
790 GMX_THROW(InternalError("invalid -g2 value"));
792 angle = gmx_angle(v1, v2);
795 GMX_THROW(InternalError("invalid -g1 value"));
797 dh.setPoint(n, angle * RAD2DEG, bPresent);
805 Angle::finishAnalysis(int /*nframes*/)
807 AbstractAverageHistogram &averageHistogram = histogramModule_->averager();
808 averageHistogram.normalizeProbability();
809 averageHistogram.done();
820 const char AngleInfo::name[] = "gangle";
821 const char AngleInfo::shortDescription[] =
824 TrajectoryAnalysisModulePointer AngleInfo::create()
826 return TrajectoryAnalysisModulePointer(new Angle);
829 } // namespace analysismodules