2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2011-2018, The GROMACS development team.
5 * Copyright (c) 2019,2020,2021, by the GROMACS development team, led by
6 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
7 * and including many others, as listed in the AUTHORS file in the
8 * top-level source directory and at http://www.gromacs.org.
10 * GROMACS is free software; you can redistribute it and/or
11 * modify it under the terms of the GNU Lesser General Public License
12 * as published by the Free Software Foundation; either version 2.1
13 * of the License, or (at your option) any later version.
15 * GROMACS is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 * Lesser General Public License for more details.
20 * You should have received a copy of the GNU Lesser General Public
21 * License along with GROMACS; if not, see
22 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
23 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
25 * If you want to redistribute modifications to GROMACS, please
26 * consider that scientific software is very special. Version
27 * control is crucial - bugs must be traceable. We will be happy to
28 * consider code for inclusion in the official distribution, but
29 * derived work must not be called official GROMACS. Details are found
30 * in the README & COPYING files - if they are missing, get the
31 * official version at http://www.gromacs.org.
33 * To help us fund GROMACS development, we humbly ask that you cite
34 * the research papers on the package. Check out http://www.gromacs.org.
38 * Implements gmx::analysismodules::Angle.
40 * \author Teemu Murtola <teemu.murtola@gmail.com>
41 * \ingroup module_trajectoryanalysis
52 #include "gromacs/analysisdata/analysisdata.h"
53 #include "gromacs/analysisdata/modules/average.h"
54 #include "gromacs/analysisdata/modules/histogram.h"
55 #include "gromacs/analysisdata/modules/plot.h"
56 #include "gromacs/math/units.h"
57 #include "gromacs/math/utilities.h"
58 #include "gromacs/math/vec.h"
59 #include "gromacs/options/basicoptions.h"
60 #include "gromacs/options/filenameoption.h"
61 #include "gromacs/options/ioptionscontainer.h"
62 #include "gromacs/pbcutil/pbc.h"
63 #include "gromacs/selection/selection.h"
64 #include "gromacs/selection/selectionoption.h"
65 #include "gromacs/trajectory/trajectoryframe.h"
66 #include "gromacs/trajectoryanalysis/analysissettings.h"
67 #include "gromacs/utility/arrayref.h"
68 #include "gromacs/utility/classhelpers.h"
69 #include "gromacs/utility/enumerationhelpers.h"
70 #include "gromacs/utility/exceptions.h"
71 #include "gromacs/utility/gmxassert.h"
72 #include "gromacs/utility/stringutil.h"
77 namespace analysismodules
83 /********************************************************************
88 * Helper to encapsulate logic for looping over input selections.
90 * This class provides two-dimensional iteration:
91 * - Over _angle groups_, corresponding to an input selection. If the input
92 * selection list contains a single selection, that selection gets used
93 * for all angle groups.
94 * - Within a group, over _values_, each consisting of a fixed number of
95 * selection positions. If there is only a single value within a selection,
96 * that value is returned over and over again.
97 * This transparently provides the semantics of using a single selection/vector
98 * to compute angles against multiple selections/vectors as described in the
101 * This class isn't perferctly self-contained and requires the caller to know
102 * some of the internals to use it properly, but it serves its purpose for this
103 * single analysis tool by simplifying the loops.
104 * Some methods have also been tailored to allow the caller to use it a bit
107 * \ingroup module_trajectoryanalysis
109 class AnglePositionIterator
113 * Creates an iterator to loop over input selection positions.
115 * \param[in] selections List of selections.
116 * \param[in] posCountPerValue Number of selection positions that
117 * constitute a single value for the iteration.
119 * If \p selections is empty, and/or \p posCountPerValue is zero, the
120 * iterator can still be advanced and hasValue()/hasSingleValue()
121 * called, but values cannot be accessed.
123 AnglePositionIterator(const SelectionList& selections, int posCountPerValue) :
124 selections_(selections), posCountPerValue_(posCountPerValue), currentSelection_(0), nextPosition_(0)
128 //! Advances the iterator to the next group of angles.
131 if (selections_.size() > 1)
137 //! Advances the iterator to the next angle in the current group.
140 if (!hasSingleValue())
142 nextPosition_ += posCountPerValue_;
147 * Returns whether this iterator represents any values.
149 * If the return value is `false`, only nextGroup(), nextValue() and
150 * hasSingleValue() are allowed to be called.
152 bool hasValue() const { return !selections_.empty(); }
154 * Returns whether the current selection only contains a single value.
156 * Returns `false` if hasValue() returns false, which allows cutting
157 * some corners in consistency checks.
159 bool hasSingleValue() const
161 return hasValue() && currentSelection().posCount() == posCountPerValue_;
163 //! Returns whether the current selection is dynamic.
164 bool isDynamic() const { 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& currentSelection() const
199 GMX_ASSERT(currentSelection_ < ssize(selections_), "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, "Accessing positions for an invalid angle type");
218 GMX_ASSERT(nextPosition_ + posCountPerValue_ <= currentSelection().posCount(),
219 "Accessing an invalid position");
220 for (int i = 0; i < posCountPerValue_; ++i)
222 copy_rvec(currentPosition(i).x(), x[i]);
227 const SelectionList& selections_;
228 const int posCountPerValue_;
229 int currentSelection_;
232 GMX_DISALLOW_COPY_AND_ASSIGN(AnglePositionIterator);
235 /********************************************************************
236 * Actual analysis module
239 //! How to interpret the selections in -group1.
240 enum class Group1Type : int
248 //! How to interpret the selections in -group2.
249 enum class Group2Type : int
259 //! String values corresponding to Group1Type.
260 const EnumerationArray<Group1Type, const char*> c_group1TypeEnumNames = {
261 { "angle", "dihedral", "vector", "plane" }
263 //! String values corresponding to Group2Type.
264 const EnumerationArray<Group2Type, const char*> c_group2TypeEnumNames = {
265 { "none", "vector", "plane", "t0", "z", "sphnorm" }
268 class Angle : public TrajectoryAnalysisModule
273 void initOptions(IOptionsContainer* options, TrajectoryAnalysisSettings* settings) override;
274 void optionsFinished(TrajectoryAnalysisSettings* settings) override;
275 void initAnalysis(const TrajectoryAnalysisSettings& settings, const TopologyInformation& top) override;
277 void analyzeFrame(int frnr, const t_trxframe& fr, t_pbc* pbc, TrajectoryAnalysisModuleData* pdata) override;
279 void finishAnalysis(int nframes) override;
280 void writeOutput() override;
283 void initFromSelections(const SelectionList& sel1, const SelectionList& sel2);
284 void checkSelections(const SelectionList& sel1, const SelectionList& sel2) const;
288 SelectionOptionInfo* sel1info_;
289 SelectionOptionInfo* sel2info_;
290 std::string fnAverage_;
292 std::string fnHistogram_;
298 AnalysisData angles_;
299 AnalysisDataFrameAverageModulePointer averageModule_;
300 AnalysisDataSimpleHistogramModulePointer histogramModule_;
302 std::vector<int> angleCount_;
305 std::vector<std::vector<RVec>> vt0_;
307 // Copy and assign disallowed by base.
313 g1type_(Group1Type::Angle),
314 g2type_(Group2Type::None),
319 averageModule_ = std::make_unique<AnalysisDataFrameAverageModule>();
320 angles_.addModule(averageModule_);
321 histogramModule_ = std::make_unique<AnalysisDataSimpleHistogramModule>();
322 angles_.addModule(histogramModule_);
324 registerAnalysisDataset(&angles_, "angle");
325 registerBasicDataset(averageModule_.get(), "average");
326 registerBasicDataset(&histogramModule_->averager(), "histogram");
330 void Angle::initOptions(IOptionsContainer* options, TrajectoryAnalysisSettings* settings)
332 static const char* const desc[] = {
333 "[THISMODULE] computes different types of angles between vectors.",
334 "It supports both vectors defined by two positions and normals of",
335 "planes defined by three positions.",
336 "The z axis or the local normal of a sphere can also be used as",
337 "one of the vectors.",
338 "There are also convenience options 'angle' and 'dihedral' for",
339 "calculating bond angles and dihedrals defined by three/four",
341 "The type of the angle is specified with [TT]-g1[tt] and [TT]-g2[tt].",
342 "If [TT]-g1[tt] is [TT]angle[tt] or [TT]dihedral[tt], [TT]-g2[tt]",
343 "should not be specified.",
344 "In this case, [TT]-group1[tt] should specify one or more selections,",
345 "and each should contain triplets or quartets of positions that define",
346 "the angles to be calculated.[PAR]",
347 "If [TT]-g1[tt] is [TT]vector[tt] or [TT]plane[tt], [TT]-group1[tt]",
348 "should specify selections that contain either pairs ([TT]vector[tt])",
349 "or triplets ([TT]plane[tt]) of positions. For vectors, the positions",
350 "set the endpoints of the vector, and for planes, the three positions",
351 "are used to calculate the normal of the plane. In both cases,",
352 "[TT]-g2[tt] specifies the other vector to use (see below).[PAR]",
353 "With [TT]-g2 vector[tt] or [TT]-g2 plane[tt], [TT]-group2[tt] should",
354 "specify another set of vectors. [TT]-group1[tt] and [TT]-group2[tt]",
355 "should specify the same number of selections. It is also allowed to",
356 "only have a single selection for one of the options, in which case",
357 "the same selection is used with each selection in the other group.",
358 "Similarly, for each selection in [TT]-group1[tt], the corresponding",
359 "selection in [TT]-group2[tt] should specify the same number of",
360 "vectors or a single vector. In the latter case, the angle is",
361 "calculated between that single vector and each vector from the other",
363 "With [TT]-g2 sphnorm[tt], each selection in [TT]-group2[tt] should",
364 "specify a single position that is the center of the sphere.",
365 "The second vector is calculated as the vector from the center to the",
366 "midpoint of the positions specified by [TT]-group1[tt].[PAR]",
367 "With [TT]-g2 z[tt], [TT]-group2[tt] is not necessary, and angles",
368 "between the first vectors and the positive Z axis are calculated.[PAR]",
369 "With [TT]-g2 t0[tt], [TT]-group2[tt] is not necessary, and angles",
370 "are calculated from the vectors as they are in the first frame.[PAR]",
371 "There are three options for output:",
372 "[TT]-oav[tt] writes an xvg file with the time and the average angle",
374 "[TT]-oall[tt] writes all the individual angles.",
375 "[TT]-oh[tt] writes a histogram of the angles. The bin width can be",
376 "set with [TT]-binw[tt].",
377 "For [TT]-oav[tt] and [TT]-oh[tt], separate average/histogram is",
378 "computed for each selection in [TT]-group1[tt]."
381 settings->setHelpText(desc);
383 options->addOption(FileNameOption("oav")
384 .filetype(OptionFileType::Plot)
387 .defaultBasename("angaver")
388 .description("Average angles as a function of time"));
389 options->addOption(FileNameOption("oall")
390 .filetype(OptionFileType::Plot)
393 .defaultBasename("angles")
394 .description("All angles as a function of time"));
395 options->addOption(FileNameOption("oh")
396 .filetype(OptionFileType::Plot)
398 .store(&fnHistogram_)
399 .defaultBasename("anghist")
400 .description("Histogram of the angles"));
402 options->addOption(EnumOption<Group1Type>("g1")
403 .enumValue(c_group1TypeEnumNames)
405 .description("Type of analysis/first vector group"));
407 EnumOption<Group2Type>("g2").enumValue(c_group2TypeEnumNames).store(&g2type_).description("Type of second vector group"));
409 DoubleOption("binw").store(&binWidth_).description("Binwidth for -oh in degrees"));
411 sel1info_ = options->addOption(
412 SelectionOption("group1").required().dynamicMask().storeVector(&sel1_).multiValue().description(
413 "First analysis/vector selection"));
414 sel2info_ = options->addOption(
415 SelectionOption("group2").dynamicMask().storeVector(&sel2_).multiValue().description(
416 "Second analysis/vector selection"));
420 void Angle::optionsFinished(TrajectoryAnalysisSettings* /* settings */)
422 const bool bSingle = (g1type_ == Group1Type::Angle || g1type_ == Group1Type::Dihedral);
424 if (bSingle && g2type_ != Group2Type::None)
427 InconsistentInputError("Cannot use a second group (-g2) with "
428 "-g1 angle or dihedral"));
430 if (bSingle && sel2info_->isSet())
433 InconsistentInputError("Cannot provide a second selection "
434 "(-group2) with -g1 angle or dihedral"));
436 if (!bSingle && g2type_ == Group2Type::None)
439 InconsistentInputError("Should specify a second group (-g2) "
440 "if the first group is not an angle or a dihedral"));
443 // Set up the number of positions per angle.
446 case Group1Type::Angle: natoms1_ = 3; break;
447 case Group1Type::Dihedral: natoms1_ = 4; break;
448 case Group1Type::Vector: natoms1_ = 2; break;
449 case Group1Type::Plane: natoms1_ = 3; break;
450 default: GMX_THROW(InternalError("invalid -g1 value"));
454 case Group2Type::None: natoms2_ = 0; break;
455 case Group2Type::Vector: natoms2_ = 2; break;
456 case Group2Type::Plane: natoms2_ = 3; break;
457 case Group2Type::TimeZero: // Intended to fall through
458 case Group2Type::Z: natoms2_ = 0; break;
459 case Group2Type::SphereNormal: natoms2_ = 1; break;
460 default: GMX_THROW(InternalError("invalid -g2 value"));
462 if (natoms2_ == 0 && sel2info_->isSet())
464 GMX_THROW(InconsistentInputError(
465 "Cannot provide a second selection (-group2) with -g2 t0 or z"));
467 // TODO: If bSingle is not set, the second selection option should be
472 void Angle::initFromSelections(const SelectionList& sel1, const SelectionList& sel2)
474 const int angleGroups = std::max(sel1.size(), sel2.size());
475 const bool bHasSecondSelection = natoms2_ > 0;
477 if (bHasSecondSelection && sel1.size() != sel2.size() && std::min(sel1.size(), sel2.size()) != 1)
479 GMX_THROW(InconsistentInputError(
480 "-group1 and -group2 should specify the same number of selections"));
483 AnglePositionIterator iter1(sel1, natoms1_);
484 AnglePositionIterator iter2(sel2, natoms2_);
485 for (int g = 0; g < angleGroups; ++g, iter1.nextGroup(), iter2.nextGroup())
487 const int posCount1 = iter1.currentSelection().posCount();
488 if (natoms1_ > 1 && posCount1 % natoms1_ != 0)
490 GMX_THROW(InconsistentInputError(formatString(
491 "Number of positions in selection %d in the first group not divisible by %d",
492 static_cast<int>(g + 1),
495 const int angleCount1 = posCount1 / natoms1_;
496 int angleCount = angleCount1;
498 if (bHasSecondSelection)
500 const int posCount2 = iter2.currentSelection().posCount();
501 if (natoms2_ > 1 && posCount2 % natoms2_ != 0)
503 GMX_THROW(InconsistentInputError(
504 formatString("Number of positions in selection %d in the second group not "
506 static_cast<int>(g + 1),
509 if (g2type_ == Group2Type::SphereNormal && posCount2 != 1)
511 GMX_THROW(InconsistentInputError(
512 "The second group should contain a single position with -g2 sphnorm"));
515 const int angleCount2 = posCount2 / natoms2_;
516 angleCount = std::max(angleCount1, angleCount2);
517 if (angleCount1 != angleCount2 && std::min(angleCount1, angleCount2) != 1)
519 GMX_THROW(InconsistentInputError(
520 "Number of vectors defined by the two groups are not the same"));
523 angleCount_.push_back(angleCount);
528 void Angle::checkSelections(const SelectionList& sel1, const SelectionList& sel2) const
530 AnglePositionIterator iter1(sel1, natoms1_);
531 AnglePositionIterator iter2(sel2, natoms2_);
532 for (size_t g = 0; g < angleCount_.size(); ++g, iter1.nextGroup(), iter2.nextGroup())
534 if (iter1.isDynamic() || (iter2.hasValue() && iter2.isDynamic()))
536 for (int n = 0; n < angleCount_[g]; ++n, iter1.nextValue(), iter2.nextValue())
539 if (!iter1.allValuesConsistentlySelected())
543 if (!iter2.allValuesConsistentlySelected())
547 if (angleCount_[g] > 1)
549 if (iter1.hasSingleValue() && !iter1.currentValuesSelected())
553 if (iter2.hasSingleValue() && !iter2.currentValuesSelected())
559 && (angleCount_[g] == 1 || (!iter1.hasSingleValue() && !iter2.hasSingleValue()))
560 && iter1.currentValuesSelected() != iter2.currentValuesSelected())
566 std::string message = formatString(
567 "Dynamic selection %d does not select "
568 "a consistent set of angles over the frames",
569 static_cast<int>(g + 1));
570 GMX_THROW(InconsistentInputError(message));
578 void Angle::initAnalysis(const TrajectoryAnalysisSettings& settings, const TopologyInformation& /* top */)
580 initFromSelections(sel1_, sel2_);
582 // checkSelections() ensures that both selection lists have the same size.
583 angles_.setDataSetCount(angleCount_.size());
584 for (size_t i = 0; i < angleCount_.size(); ++i)
586 angles_.setColumnCount(i, angleCount_[i]);
588 double histogramMin = (g1type_ == Group1Type::Dihedral ? -180.0 : 0);
589 histogramModule_->init(histogramFromRange(histogramMin, 180.0).binWidth(binWidth_).includeAll());
591 if (g2type_ == Group2Type::TimeZero)
593 vt0_.resize(sel1_.size());
594 for (size_t g = 0; g < sel1_.size(); ++g)
596 vt0_[g].resize(sel1_[g].posCount() / natoms1_);
600 if (!fnAverage_.empty())
602 AnalysisDataPlotModulePointer plotm(new AnalysisDataPlotModule(settings.plotSettings()));
603 plotm->setFileName(fnAverage_);
604 plotm->setTitle("Average angle");
605 plotm->setXAxisIsTime();
606 plotm->setYLabel("Angle (degrees)");
607 // TODO: Consider adding information about the second selection,
608 // and/or a subtitle describing what kind of angle this is.
609 for (size_t g = 0; g < sel1_.size(); ++g)
611 plotm->appendLegend(sel1_[g].name());
613 averageModule_->addModule(plotm);
618 AnalysisDataPlotModulePointer plotm(new AnalysisDataPlotModule(settings.plotSettings()));
619 plotm->setFileName(fnAll_);
620 plotm->setTitle("Angle");
621 plotm->setXAxisIsTime();
622 plotm->setYLabel("Angle (degrees)");
623 // TODO: Add legends? (there can be a massive amount of columns)
624 angles_.addModule(plotm);
627 if (!fnHistogram_.empty())
629 AnalysisDataPlotModulePointer plotm(new AnalysisDataPlotModule(settings.plotSettings()));
630 plotm->setFileName(fnHistogram_);
631 plotm->setTitle("Angle histogram");
632 plotm->setXLabel("Angle (degrees)");
633 plotm->setYLabel("Probability");
634 // TODO: Consider adding information about the second selection,
635 // and/or a subtitle describing what kind of angle this is.
636 for (size_t g = 0; g < sel1_.size(); ++g)
638 plotm->appendLegend(sel1_[g].name());
640 histogramModule_->averager().addModule(plotm);
645 //! Helper method to calculate a vector from two or three positions.
646 void calc_vec(int natoms, rvec x[], t_pbc* pbc, rvec xout, rvec cout)
653 pbc_dx(pbc, x[1], x[0], xout);
657 rvec_sub(x[1], x[0], xout);
659 svmul(0.5, xout, cout);
660 rvec_add(x[0], cout, cout);
667 pbc_dx(pbc, x[1], x[0], v1);
668 pbc_dx(pbc, x[2], x[0], v2);
672 rvec_sub(x[1], x[0], v1);
673 rvec_sub(x[2], x[0], v2);
676 rvec_add(x[0], x[1], cout);
677 rvec_add(cout, x[2], cout);
678 svmul(1.0 / 3.0, cout, cout);
681 default: GMX_RELEASE_ASSERT(false, "Incorrectly initialized number of atoms");
686 void Angle::analyzeFrame(int frnr, const t_trxframe& fr, t_pbc* pbc, TrajectoryAnalysisModuleData* pdata)
688 AnalysisDataHandle dh = pdata->dataHandle(angles_);
689 const SelectionList& sel1 = pdata->parallelSelections(sel1_);
690 const SelectionList& sel2 = pdata->parallelSelections(sel2_);
692 checkSelections(sel1, sel2);
694 dh.startFrame(frnr, fr.time);
696 AnglePositionIterator iter1(sel1, natoms1_);
697 AnglePositionIterator iter2(sel2, natoms2_);
698 for (size_t g = 0; g < angleCount_.size(); ++g, iter1.nextGroup(), iter2.nextGroup())
703 // v2 & c2 are conditionally set in the switch statement below, and conditionally
704 // used in a different switch statement later. Apparently the clang static analyzer
705 // thinks there are cases where they can be used uninitialzed (which I cannot find),
706 // but to avoid trouble if we ever change just one of the switch statements it
707 // makes sense to clear them outside the first switch.
714 case Group2Type::Z: v2[ZZ] = 1.0; break;
715 case Group2Type::SphereNormal: copy_rvec(sel2_[g].position(0).x(), c2); break;
722 for (int n = 0; n < angleCount_[g]; ++n, iter1.nextValue(), iter2.nextValue())
725 // x[] will be assigned below based on the number of atoms used to initialize iter1,
726 // which in turn should correspond perfectly to g1type_ (which determines how many we
727 // read), but unsurprisingly the static analyzer chokes a bit on that.
731 // checkSelections() ensures that this reflects all the involved
733 const bool bPresent = iter1.currentValuesSelected() && iter2.currentValuesSelected();
734 iter1.getCurrentPositions(x);
737 case Group1Type::Angle:
740 pbc_dx(pbc, x[0], x[1], v1);
741 pbc_dx(pbc, x[2], x[1], v2);
745 rvec_sub(x[0], x[1], v1);
746 rvec_sub(x[2], x[1], v2);
748 angle = gmx_angle(v1, v2);
750 case Group1Type::Dihedral:
755 pbc_dx(pbc, x[0], x[1], dx[0]);
756 pbc_dx(pbc, x[2], x[1], dx[1]);
757 pbc_dx(pbc, x[2], x[3], dx[2]);
761 rvec_sub(x[0], x[1], dx[0]);
762 rvec_sub(x[2], x[1], dx[1]);
763 rvec_sub(x[2], x[3], dx[2]);
765 cprod(dx[0], dx[1], v1);
766 cprod(dx[1], dx[2], v2);
767 angle = gmx_angle(v1, v2);
768 real ipr = iprod(dx[0], v2);
775 case Group1Type::Vector:
776 case Group1Type::Plane:
777 calc_vec(natoms1_, x, pbc, v1, c1);
780 case Group2Type::Vector:
781 case Group2Type::Plane:
782 iter2.getCurrentPositions(x);
783 calc_vec(natoms2_, x, pbc, v2, c2);
785 case Group2Type::TimeZero:
786 // FIXME: This is not parallelizable.
789 copy_rvec(v1, vt0_[g][n]);
791 copy_rvec(vt0_[g][n], v2);
793 case Group2Type::Z: c1[XX] = c1[YY] = 0.0; break;
794 case Group2Type::SphereNormal:
797 pbc_dx(pbc, c1, c2, v2);
801 rvec_sub(c1, c2, v2);
804 default: GMX_THROW(InternalError("invalid -g2 value"));
806 angle = gmx_angle(v1, v2);
808 default: GMX_THROW(InternalError("invalid -g1 value"));
810 dh.setPoint(n, angle * gmx::c_rad2Deg, bPresent);
817 void Angle::finishAnalysis(int /*nframes*/)
819 AbstractAverageHistogram& averageHistogram = histogramModule_->averager();
820 averageHistogram.normalizeProbability();
821 averageHistogram.done();
825 void Angle::writeOutput() {}
829 const char AngleInfo::name[] = "gangle";
830 const char AngleInfo::shortDescription[] = "Calculate angles";
832 TrajectoryAnalysisModulePointer AngleInfo::create()
834 return TrajectoryAnalysisModulePointer(new Angle);
837 } // namespace analysismodules