Tidy: modernize-use-nullptr
[gromacs.git] / src / gromacs / trajectoryanalysis / modules / angle.cpp
blob058e5329c617a98ce8ada6459f5399acf6b43ee6
1 /*
2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2011,2012,2013,2014,2015,2016,2017, 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.
35 /*! \internal \file
36 * \brief
37 * Implements gmx::analysismodules::Angle.
39 * \author Teemu Murtola <teemu.murtola@gmail.com>
40 * \ingroup module_trajectoryanalysis
42 #include "gmxpre.h"
44 #include "angle.h"
46 #include <algorithm>
47 #include <string>
48 #include <vector>
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/math/units.h"
55 #include "gromacs/math/vec.h"
56 #include "gromacs/options/basicoptions.h"
57 #include "gromacs/options/filenameoption.h"
58 #include "gromacs/options/ioptionscontainer.h"
59 #include "gromacs/pbcutil/pbc.h"
60 #include "gromacs/selection/selection.h"
61 #include "gromacs/selection/selectionoption.h"
62 #include "gromacs/trajectory/trajectoryframe.h"
63 #include "gromacs/trajectoryanalysis/analysissettings.h"
64 #include "gromacs/utility/arrayref.h"
65 #include "gromacs/utility/exceptions.h"
66 #include "gromacs/utility/gmxassert.h"
67 #include "gromacs/utility/stringutil.h"
69 namespace gmx
72 namespace analysismodules
75 namespace
78 /********************************************************************
79 * Helper classes
82 /*! \brief
83 * Helper to encapsulate logic for looping over input selections.
85 * This class provides two-dimensional iteration:
86 * - Over _angle groups_, corresponding to an input selection. If the input
87 * selection list contains a single selection, that selection gets used
88 * for all angle groups.
89 * - Within a group, over _values_, each consisting of a fixed number of
90 * selection positions. If there is only a single value within a selection,
91 * that value is returned over and over again.
92 * This transparently provides the semantics of using a single selection/vector
93 * to compute angles against multiple selections/vectors as described in the
94 * tool help text.
96 * This class isn't perferctly self-contained and requires the caller to know
97 * some of the internals to use it properly, but it serves its purpose for this
98 * single analysis tool by simplifying the loops.
99 * Some methods have also been tailored to allow the caller to use it a bit
100 * more easily.
102 * \ingroup module_trajectoryanalysis
104 class AnglePositionIterator
106 public:
107 /*! \brief
108 * Creates an iterator to loop over input selection positions.
110 * \param[in] selections List of selections.
111 * \param[in] posCountPerValue Number of selection positions that
112 * constitute a single value for the iteration.
114 * If \p selections is empty, and/or \p posCountPerValue is zero, the
115 * iterator can still be advanced and hasValue()/hasSingleValue()
116 * called, but values cannot be accessed.
118 AnglePositionIterator(const SelectionList &selections,
119 int posCountPerValue)
120 : selections_(selections), posCountPerValue_(posCountPerValue),
121 currentSelection_(0), nextPosition_(0)
125 //! Advances the iterator to the next group of angles.
126 void nextGroup()
128 if (selections_.size() > 1)
130 ++currentSelection_;
132 nextPosition_ = 0;
134 //! Advances the iterator to the next angle in the current group.
135 void nextValue()
137 if (!hasSingleValue())
139 nextPosition_ += posCountPerValue_;
143 /*! \brief
144 * Returns whether this iterator represents any values.
146 * If the return value is `false`, only nextGroup(), nextValue() and
147 * hasSingleValue() are allowed to be called.
149 bool hasValue() const
151 return !selections_.empty();
153 /*! \brief
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
166 return currentSelection().isDynamic();
168 /*! \brief
169 * Returns whether positions in the current value are either all
170 * selected or all unselected.
172 bool allValuesConsistentlySelected() const
174 if (posCountPerValue_ <= 1)
176 return true;
178 const bool bSelected = currentPosition(0).selected();
179 for (int i = 1; i < posCountPerValue_; ++i)
181 if (currentPosition(i).selected() != bSelected)
183 return false;
186 return true;
188 /*! \brief
189 * Returns whether positions in the current value are selected.
191 * Only works reliably if allValuesConsistentlySelected() returns
192 * `true`.
194 bool currentValuesSelected() const
196 return selections_.empty() || currentPosition(0).selected();
199 //! Returns the currently active selection.
200 const Selection &currentSelection() const
202 GMX_ASSERT(currentSelection_ < static_cast<int>(selections_.size()),
203 "Accessing an invalid selection");
204 return selections_[currentSelection_];
206 //! Returns the `i`th position for the current value.
207 SelectionPosition currentPosition(int i) const
209 return currentSelection().position(nextPosition_ + i);
211 /*! \brief
212 * Extracts all coordinates corresponding to the current value.
214 * \param[out] x Array to which the positions are extracted.
216 * \p x should contain at minimum the number of positions per value
217 * passed to the constructor.
219 void getCurrentPositions(rvec x[]) const
221 GMX_ASSERT(posCountPerValue_ > 0,
222 "Accessing positions for an invalid angle type");
223 GMX_ASSERT(nextPosition_ + posCountPerValue_ <= currentSelection().posCount(),
224 "Accessing an invalid position");
225 for (int i = 0; i < posCountPerValue_; ++i)
227 copy_rvec(currentPosition(i).x(), x[i]);
231 private:
232 const SelectionList &selections_;
233 const int posCountPerValue_;
234 int currentSelection_;
235 int nextPosition_;
237 GMX_DISALLOW_COPY_AND_ASSIGN(AnglePositionIterator);
240 /********************************************************************
241 * Actual analysis module
244 //! How to interpret the selections in -group1.
245 enum Group1Type
247 Group1Type_Angle,
248 Group1Type_Dihedral,
249 Group1Type_Vector,
250 Group1Type_Plane
252 //! How to interpret the selections in -group2.
253 enum Group2Type
255 Group2Type_None,
256 Group2Type_Vector,
257 Group2Type_Plane,
258 Group2Type_TimeZero,
259 Group2Type_Z,
260 Group2Type_SphereNormal
262 //! String values corresponding to Group1Type.
263 const char *const cGroup1TypeEnum[] =
264 { "angle", "dihedral", "vector", "plane" };
265 //! String values corresponding to Group2Type.
266 const char *const cGroup2TypeEnum[] =
267 { "none", "vector", "plane", "t0", "z", "sphnorm" };
269 class Angle : public TrajectoryAnalysisModule
271 public:
272 Angle();
274 virtual void initOptions(IOptionsContainer *options,
275 TrajectoryAnalysisSettings *settings);
276 virtual void optionsFinished(TrajectoryAnalysisSettings *settings);
277 virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
278 const TopologyInformation &top);
280 virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
281 TrajectoryAnalysisModuleData *pdata);
283 virtual void finishAnalysis(int nframes);
284 virtual void writeOutput();
286 private:
287 void initFromSelections(const SelectionList &sel1,
288 const SelectionList &sel2);
289 void checkSelections(const SelectionList &sel1,
290 const SelectionList &sel2) const;
292 SelectionList sel1_;
293 SelectionList sel2_;
294 SelectionOptionInfo *sel1info_;
295 SelectionOptionInfo *sel2info_;
296 std::string fnAverage_;
297 std::string fnAll_;
298 std::string fnHistogram_;
300 Group1Type g1type_;
301 Group2Type g2type_;
302 double binWidth_;
304 AnalysisData angles_;
305 AnalysisDataFrameAverageModulePointer averageModule_;
306 AnalysisDataSimpleHistogramModulePointer histogramModule_;
308 std::vector<int> angleCount_;
309 int natoms1_;
310 int natoms2_;
311 std::vector<std::vector<RVec> > vt0_;
313 // Copy and assign disallowed by base.
316 Angle::Angle()
317 : sel1info_(nullptr), sel2info_(nullptr),
318 g1type_(Group1Type_Angle), g2type_(Group2Type_None),
319 binWidth_(1.0), natoms1_(0), natoms2_(0)
321 averageModule_.reset(new AnalysisDataFrameAverageModule());
322 angles_.addModule(averageModule_);
323 histogramModule_.reset(new AnalysisDataSimpleHistogramModule());
324 angles_.addModule(histogramModule_);
326 registerAnalysisDataset(&angles_, "angle");
327 registerBasicDataset(averageModule_.get(), "average");
328 registerBasicDataset(&histogramModule_->averager(), "histogram");
332 void
333 Angle::initOptions(IOptionsContainer *options, TrajectoryAnalysisSettings *settings)
335 static const char *const desc[] = {
336 "[THISMODULE] computes different types of angles between vectors.",
337 "It supports both vectors defined by two positions and normals of",
338 "planes defined by three positions.",
339 "The z axis or the local normal of a sphere can also be used as",
340 "one of the vectors.",
341 "There are also convenience options 'angle' and 'dihedral' for",
342 "calculating bond angles and dihedrals defined by three/four",
343 "positions.[PAR]",
344 "The type of the angle is specified with [TT]-g1[tt] and [TT]-g2[tt].",
345 "If [TT]-g1[tt] is [TT]angle[tt] or [TT]dihedral[tt], [TT]-g2[tt]",
346 "should not be specified.",
347 "In this case, [TT]-group1[tt] should specify one or more selections,",
348 "and each should contain triplets or quartets of positions that define",
349 "the angles to be calculated.[PAR]",
350 "If [TT]-g1[tt] is [TT]vector[tt] or [TT]plane[tt], [TT]-group1[tt]",
351 "should specify selections that contain either pairs ([TT]vector[tt])",
352 "or triplets ([TT]plane[tt]) of positions. For vectors, the positions",
353 "set the endpoints of the vector, and for planes, the three positions",
354 "are used to calculate the normal of the plane. In both cases,",
355 "[TT]-g2[tt] specifies the other vector to use (see below).[PAR]",
356 "With [TT]-g2 vector[tt] or [TT]-g2 plane[tt], [TT]-group2[tt] should",
357 "specify another set of vectors. [TT]-group1[tt] and [TT]-group2[tt]",
358 "should specify the same number of selections. It is also allowed to",
359 "only have a single selection for one of the options, in which case",
360 "the same selection is used with each selection in the other group.",
361 "Similarly, for each selection in [TT]-group1[tt], the corresponding",
362 "selection in [TT]-group2[tt] should specify the same number of",
363 "vectors or a single vector. In the latter case, the angle is",
364 "calculated between that single vector and each vector from the other",
365 "selection.[PAR]",
366 "With [TT]-g2 sphnorm[tt], each selection in [TT]-group2[tt] should",
367 "specify a single position that is the center of the sphere.",
368 "The second vector is calculated as the vector from the center to the",
369 "midpoint of the positions specified by [TT]-group1[tt].[PAR]",
370 "With [TT]-g2 z[tt], [TT]-group2[tt] is not necessary, and angles",
371 "between the first vectors and the positive Z axis are calculated.[PAR]",
372 "With [TT]-g2 t0[tt], [TT]-group2[tt] is not necessary, and angles",
373 "are calculated from the vectors as they are in the first frame.[PAR]",
374 "There are three options for output:",
375 "[TT]-oav[tt] writes an xvg file with the time and the average angle",
376 "for each frame.",
377 "[TT]-oall[tt] writes all the individual angles.",
378 "[TT]-oh[tt] writes a histogram of the angles. The bin width can be",
379 "set with [TT]-binw[tt].",
380 "For [TT]-oav[tt] and [TT]-oh[tt], separate average/histogram is",
381 "computed for each selection in [TT]-group1[tt]."
384 settings->setHelpText(desc);
386 options->addOption(FileNameOption("oav").filetype(eftPlot).outputFile()
387 .store(&fnAverage_).defaultBasename("angaver")
388 .description("Average angles as a function of time"));
389 options->addOption(FileNameOption("oall").filetype(eftPlot).outputFile()
390 .store(&fnAll_).defaultBasename("angles")
391 .description("All angles as a function of time"));
392 options->addOption(FileNameOption("oh").filetype(eftPlot).outputFile()
393 .store(&fnHistogram_).defaultBasename("anghist")
394 .description("Histogram of the angles"));
396 options->addOption(EnumOption<Group1Type>("g1").enumValue(cGroup1TypeEnum)
397 .store(&g1type_)
398 .description("Type of analysis/first vector group"));
399 options->addOption(EnumOption<Group2Type>("g2").enumValue(cGroup2TypeEnum)
400 .store(&g2type_)
401 .description("Type of second vector group"));
402 options->addOption(DoubleOption("binw").store(&binWidth_)
403 .description("Binwidth for -oh in degrees"));
405 sel1info_ = options->addOption(SelectionOption("group1")
406 .required().dynamicMask().storeVector(&sel1_)
407 .multiValue()
408 .description("First analysis/vector selection"));
409 sel2info_ = options->addOption(SelectionOption("group2")
410 .dynamicMask().storeVector(&sel2_)
411 .multiValue()
412 .description("Second analysis/vector selection"));
416 void
417 Angle::optionsFinished(TrajectoryAnalysisSettings * /* settings */)
419 const bool bSingle = (g1type_ == Group1Type_Angle || g1type_ == Group1Type_Dihedral);
421 if (bSingle && g2type_ != Group2Type_None)
423 GMX_THROW(InconsistentInputError("Cannot use a second group (-g2) with "
424 "-g1 angle or dihedral"));
426 if (bSingle && sel2info_->isSet())
428 GMX_THROW(InconsistentInputError("Cannot provide a second selection "
429 "(-group2) with -g1 angle or dihedral"));
431 if (!bSingle && g2type_ == Group2Type_None)
433 GMX_THROW(InconsistentInputError("Should specify a second group (-g2) "
434 "if the first group is not an angle or a dihedral"));
437 // Set up the number of positions per angle.
438 switch (g1type_)
440 case Group1Type_Angle: natoms1_ = 3; break;
441 case Group1Type_Dihedral: natoms1_ = 4; break;
442 case Group1Type_Vector: natoms1_ = 2; break;
443 case Group1Type_Plane: natoms1_ = 3; break;
444 default:
445 GMX_THROW(InternalError("invalid -g1 value"));
447 switch (g2type_)
449 case Group2Type_None: natoms2_ = 0; break;
450 case Group2Type_Vector: natoms2_ = 2; break;
451 case Group2Type_Plane: natoms2_ = 3; break;
452 case Group2Type_TimeZero: natoms2_ = 0; break;
453 case Group2Type_Z: natoms2_ = 0; break;
454 case Group2Type_SphereNormal: natoms2_ = 1; break;
455 default:
456 GMX_THROW(InternalError("invalid -g2 value"));
458 if (natoms2_ == 0 && sel2info_->isSet())
460 GMX_THROW(InconsistentInputError("Cannot provide a second selection (-group2) with -g2 t0 or z"));
462 // TODO: If bSingle is not set, the second selection option should be
463 // required.
467 void
468 Angle::initFromSelections(const SelectionList &sel1,
469 const SelectionList &sel2)
471 const int angleGroups = std::max(sel1.size(), sel2.size());
472 const bool bHasSecondSelection = natoms2_ > 0;
474 if (bHasSecondSelection && sel1.size() != sel2.size()
475 && std::min(sel1.size(), sel2.size()) != 1)
477 GMX_THROW(InconsistentInputError(
478 "-group1 and -group2 should specify the same number of selections"));
481 AnglePositionIterator iter1(sel1, natoms1_);
482 AnglePositionIterator iter2(sel2, natoms2_);
483 for (int g = 0; g < angleGroups; ++g, iter1.nextGroup(), iter2.nextGroup())
485 const int posCount1 = iter1.currentSelection().posCount();
486 if (natoms1_ > 1 && posCount1 % natoms1_ != 0)
488 GMX_THROW(InconsistentInputError(formatString(
489 "Number of positions in selection %d in the first group not divisible by %d",
490 static_cast<int>(g + 1), natoms1_)));
492 const int angleCount1 = posCount1 / natoms1_;
493 int angleCount = angleCount1;
495 if (bHasSecondSelection)
497 const int posCount2 = iter2.currentSelection().posCount();
498 if (natoms2_ > 1 && posCount2 % natoms2_ != 0)
500 GMX_THROW(InconsistentInputError(formatString(
501 "Number of positions in selection %d in the second group not divisible by %d",
502 static_cast<int>(g + 1), natoms2_)));
504 if (g2type_ == Group2Type_SphereNormal && posCount2 != 1)
506 GMX_THROW(InconsistentInputError(
507 "The second group should contain a single position with -g2 sphnorm"));
510 const int angleCount2 = posCount2 / natoms2_;
511 angleCount = std::max(angleCount1, angleCount2);
512 if (angleCount1 != angleCount2
513 && std::min(angleCount1, angleCount2) != 1)
515 GMX_THROW(InconsistentInputError(
516 "Number of vectors defined by the two groups are not the same"));
519 angleCount_.push_back(angleCount);
524 void
525 Angle::checkSelections(const SelectionList &sel1,
526 const SelectionList &sel2) const
528 AnglePositionIterator iter1(sel1, natoms1_);
529 AnglePositionIterator iter2(sel2, natoms2_);
530 for (size_t g = 0; g < angleCount_.size(); ++g, iter1.nextGroup(), iter2.nextGroup())
532 if (iter1.isDynamic() || (iter2.hasValue() && iter2.isDynamic()))
534 for (int n = 0; n < angleCount_[g]; ++n, iter1.nextValue(), iter2.nextValue())
536 bool bOk = true;
537 if (!iter1.allValuesConsistentlySelected())
539 bOk = false;
541 if (!iter2.allValuesConsistentlySelected())
543 bOk = false;
545 if (angleCount_[g] > 1)
547 if (iter1.hasSingleValue() && !iter1.currentValuesSelected())
549 bOk = false;
551 if (iter2.hasSingleValue() && !iter2.currentValuesSelected())
553 bOk = false;
556 if (iter2.hasValue()
557 && (angleCount_[g] == 1
558 || (!iter1.hasSingleValue() && !iter2.hasSingleValue()))
559 && iter1.currentValuesSelected() != iter2.currentValuesSelected())
561 bOk = false;
563 if (!bOk)
565 std::string message =
566 formatString("Dynamic selection %d does not select "
567 "a consistent set of angles over the frames",
568 static_cast<int>(g + 1));
569 GMX_THROW(InconsistentInputError(message));
577 void
578 Angle::initAnalysis(const TrajectoryAnalysisSettings &settings,
579 const TopologyInformation & /* top */)
581 initFromSelections(sel1_, sel2_);
583 // checkSelections() ensures that both selection lists have the same size.
584 angles_.setDataSetCount(angleCount_.size());
585 for (size_t i = 0; i < angleCount_.size(); ++i)
587 angles_.setColumnCount(i, angleCount_[i]);
589 double histogramMin = (g1type_ == Group1Type_Dihedral ? -180.0 : 0);
590 histogramModule_->init(histogramFromRange(histogramMin, 180.0)
591 .binWidth(binWidth_).includeAll());
593 if (g2type_ == Group2Type_TimeZero)
595 vt0_.resize(sel1_.size());
596 for (size_t g = 0; g < sel1_.size(); ++g)
598 vt0_[g].resize(sel1_[g].posCount() / natoms1_);
602 if (!fnAverage_.empty())
604 AnalysisDataPlotModulePointer plotm(
605 new AnalysisDataPlotModule(settings.plotSettings()));
606 plotm->setFileName(fnAverage_);
607 plotm->setTitle("Average angle");
608 plotm->setXAxisIsTime();
609 plotm->setYLabel("Angle (degrees)");
610 // TODO: Consider adding information about the second selection,
611 // and/or a subtitle describing what kind of angle this is.
612 for (size_t g = 0; g < sel1_.size(); ++g)
614 plotm->appendLegend(sel1_[g].name());
616 averageModule_->addModule(plotm);
619 if (!fnAll_.empty())
621 AnalysisDataPlotModulePointer plotm(
622 new AnalysisDataPlotModule(settings.plotSettings()));
623 plotm->setFileName(fnAll_);
624 plotm->setTitle("Angle");
625 plotm->setXAxisIsTime();
626 plotm->setYLabel("Angle (degrees)");
627 // TODO: Add legends? (there can be a massive amount of columns)
628 angles_.addModule(plotm);
631 if (!fnHistogram_.empty())
633 AnalysisDataPlotModulePointer plotm(
634 new AnalysisDataPlotModule(settings.plotSettings()));
635 plotm->setFileName(fnHistogram_);
636 plotm->setTitle("Angle histogram");
637 plotm->setXLabel("Angle (degrees)");
638 plotm->setYLabel("Probability");
639 // TODO: Consider adding information about the second selection,
640 // and/or a subtitle describing what kind of angle this is.
641 for (size_t g = 0; g < sel1_.size(); ++g)
643 plotm->appendLegend(sel1_[g].name());
645 histogramModule_->averager().addModule(plotm);
650 //! Helper method to calculate a vector from two or three positions.
651 static void
652 calc_vec(int natoms, rvec x[], t_pbc *pbc, rvec xout, rvec cout)
654 switch (natoms)
656 case 2:
657 if (pbc)
659 pbc_dx(pbc, x[1], x[0], xout);
661 else
663 rvec_sub(x[1], x[0], xout);
665 svmul(0.5, xout, cout);
666 rvec_add(x[0], cout, cout);
667 break;
668 case 3:
670 rvec v1, v2;
671 if (pbc)
673 pbc_dx(pbc, x[1], x[0], v1);
674 pbc_dx(pbc, x[2], x[0], v2);
676 else
678 rvec_sub(x[1], x[0], v1);
679 rvec_sub(x[2], x[0], v2);
681 cprod(v1, v2, xout);
682 rvec_add(x[0], x[1], cout);
683 rvec_add(cout, x[2], cout);
684 svmul(1.0/3.0, cout, cout);
685 break;
687 default:
688 GMX_RELEASE_ASSERT(false, "Incorrectly initialized number of atoms");
693 void
694 Angle::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
695 TrajectoryAnalysisModuleData *pdata)
697 AnalysisDataHandle dh = pdata->dataHandle(angles_);
698 const SelectionList &sel1 = pdata->parallelSelections(sel1_);
699 const SelectionList &sel2 = pdata->parallelSelections(sel2_);
701 checkSelections(sel1, sel2);
703 dh.startFrame(frnr, fr.time);
705 AnglePositionIterator iter1(sel1, natoms1_);
706 AnglePositionIterator iter2(sel2, natoms2_);
707 for (size_t g = 0; g < angleCount_.size(); ++g, iter1.nextGroup(), iter2.nextGroup())
709 rvec v1, v2;
710 rvec c1, c2;
712 // v2 & c2 are conditionally set in the switch statement below, and conditionally
713 // used in a different switch statement later. Apparently the clang static analyzer
714 // thinks there are cases where they can be used uninitialzed (which I cannot find),
715 // but to avoid trouble if we ever change just one of the switch statements it
716 // makes sense to clear them outside the first switch.
718 clear_rvec(v2);
719 clear_rvec(c2);
721 switch (g2type_)
723 case Group2Type_Z:
724 v2[ZZ] = 1.0;
725 break;
726 case Group2Type_SphereNormal:
727 copy_rvec(sel2_[g].position(0).x(), c2);
728 break;
729 default:
730 // do nothing
731 break;
734 dh.selectDataSet(g);
735 for (int n = 0; n < angleCount_[g]; ++n, iter1.nextValue(), iter2.nextValue())
737 rvec x[4];
738 // x[] will be assigned below based on the number of atoms used to initialize iter1,
739 // which in turn should correspond perfectly to g1type_ (which determines how many we read),
740 // but unsurprisingly the static analyzer chokes a bit on that.
741 clear_rvecs(4, x);
743 real angle;
744 // checkSelections() ensures that this reflects all the involved
745 // positions.
746 const bool bPresent =
747 iter1.currentValuesSelected() && iter2.currentValuesSelected();
748 iter1.getCurrentPositions(x);
749 switch (g1type_)
751 case Group1Type_Angle:
752 if (pbc)
754 pbc_dx(pbc, x[0], x[1], v1);
755 pbc_dx(pbc, x[2], x[1], v2);
757 else
759 rvec_sub(x[0], x[1], v1);
760 rvec_sub(x[2], x[1], v2);
762 angle = gmx_angle(v1, v2);
763 break;
764 case Group1Type_Dihedral:
766 rvec dx[3];
767 if (pbc)
769 // cppcheck-suppress uninitvar
770 pbc_dx(pbc, x[0], x[1], dx[0]);
771 // cppcheck-suppress uninitvar
772 pbc_dx(pbc, x[2], x[1], dx[1]);
773 // cppcheck-suppress uninitvar
774 pbc_dx(pbc, x[2], x[3], dx[2]);
776 else
778 // cppcheck-suppress uninitvar
779 rvec_sub(x[0], x[1], dx[0]);
780 // cppcheck-suppress uninitvar
781 rvec_sub(x[2], x[1], dx[1]);
782 // cppcheck-suppress uninitvar
783 rvec_sub(x[2], x[3], dx[2]);
785 cprod(dx[0], dx[1], v1);
786 cprod(dx[1], dx[2], v2);
787 angle = gmx_angle(v1, v2);
788 real ipr = iprod(dx[0], v2);
789 if (ipr < 0)
791 angle = -angle;
793 break;
795 case Group1Type_Vector:
796 case Group1Type_Plane:
797 calc_vec(natoms1_, x, pbc, v1, c1);
798 switch (g2type_)
800 case Group2Type_Vector:
801 case Group2Type_Plane:
802 iter2.getCurrentPositions(x);
803 calc_vec(natoms2_, x, pbc, v2, c2);
804 break;
805 case Group2Type_TimeZero:
806 // FIXME: This is not parallelizable.
807 if (frnr == 0)
809 copy_rvec(v1, vt0_[g][n]);
811 copy_rvec(vt0_[g][n], v2);
812 break;
813 case Group2Type_Z:
814 c1[XX] = c1[YY] = 0.0;
815 break;
816 case Group2Type_SphereNormal:
817 if (pbc)
819 pbc_dx(pbc, c1, c2, v2);
821 else
823 rvec_sub(c1, c2, v2);
825 break;
826 default:
827 GMX_THROW(InternalError("invalid -g2 value"));
829 angle = gmx_angle(v1, v2);
830 break;
831 default:
832 GMX_THROW(InternalError("invalid -g1 value"));
834 dh.setPoint(n, angle * RAD2DEG, bPresent);
837 dh.finishFrame();
841 void
842 Angle::finishAnalysis(int /*nframes*/)
844 AbstractAverageHistogram &averageHistogram = histogramModule_->averager();
845 averageHistogram.normalizeProbability();
846 averageHistogram.done();
850 void
851 Angle::writeOutput()
855 } // namespace
857 const char AngleInfo::name[] = "gangle";
858 const char AngleInfo::shortDescription[] =
859 "Calculate angles";
861 TrajectoryAnalysisModulePointer AngleInfo::create()
863 return TrajectoryAnalysisModulePointer(new Angle);
866 } // namespace analysismodules
868 } // namespace gmx