src/gromacs/trajectoryanalysis/modules/angle.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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.
   9  *
  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.
  14  *
  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.
  19  *
  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.
  24  *
  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.
  32  *
  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.
  35  */
  36 /*! \internal \file
  37  * \brief
  38  * Implements gmx::analysismodules::Angle.
  39  *
  40  * \author Teemu Murtola <teemu.murtola@gmail.com>
  41  * \ingroup module_trajectoryanalysis
  42  */
  43 #include "gmxpre.h"
  44
  45 #include "angle.h"
  46
  47 #include <algorithm>
  48 #include <memory>
  49 #include <string>
  50 #include <vector>
  51
  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"
  73
  74 namespace gmx
  75 {
  76
  77 namespace analysismodules
  78 {
  79
  80 namespace
  81 {
  82
  83 /********************************************************************
  84  * Helper classes
  85  */
  86
  87 /*! \brief
  88  * Helper to encapsulate logic for looping over input selections.
  89  *
  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
  99  * tool help text.
 100  *
 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
 105  * more easily.
 106  *
 107  * \ingroup module_trajectoryanalysis
 108  */
 109 class AnglePositionIterator
 110 {
 111 public:
 112     /*! \brief
 113      * Creates an iterator to loop over input selection positions.
 114      *
 115      * \param[in] selections       List of selections.
 116      * \param[in] posCountPerValue Number of selection positions that
 117      *     constitute a single value for the iteration.
 118      *
 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.
 122      */
 123     AnglePositionIterator(const SelectionList& selections, int posCountPerValue) :
 124         selections_(selections), posCountPerValue_(posCountPerValue), currentSelection_(0), nextPosition_(0)
 125     {
 126     }
 127
 128     //! Advances the iterator to the next group of angles.
 129     void nextGroup()
 130     {
 131         if (selections_.size() > 1)
 132         {
 133             ++currentSelection_;
 134         }
 135         nextPosition_ = 0;
 136     }
 137     //! Advances the iterator to the next angle in the current group.
 138     void nextValue()
 139     {
 140         if (!hasSingleValue())
 141         {
 142             nextPosition_ += posCountPerValue_;
 143         }
 144     }
 145
 146     /*! \brief
 147      * Returns whether this iterator represents any values.
 148      *
 149      * If the return value is `false`, only nextGroup(), nextValue() and
 150      * hasSingleValue() are allowed to be called.
 151      */
 152     bool hasValue() const { return !selections_.empty(); }
 153     /*! \brief
 154      * Returns whether the current selection only contains a single value.
 155      *
 156      * Returns `false` if hasValue() returns false, which allows cutting
 157      * some corners in consistency checks.
 158      */
 159     bool hasSingleValue() const
 160     {
 161         return hasValue() && currentSelection().posCount() == posCountPerValue_;
 162     }
 163     //! Returns whether the current selection is dynamic.
 164     bool isDynamic() const { return currentSelection().isDynamic(); }
 165     /*! \brief
 166      * Returns whether positions in the current value are either all
 167      * selected or all unselected.
 168      */
 169     bool allValuesConsistentlySelected() const
 170     {
 171         if (posCountPerValue_ <= 1)
 172         {
 173             return true;
 174         }
 175         const bool bSelected = currentPosition(0).selected();
 176         for (int i = 1; i < posCountPerValue_; ++i)
 177         {
 178             if (currentPosition(i).selected() != bSelected)
 179             {
 180                 return false;
 181             }
 182         }
 183         return true;
 184     }
 185     /*! \brief
 186      * Returns whether positions in the current value are selected.
 187      *
 188      * Only works reliably if allValuesConsistentlySelected() returns
 189      * `true`.
 190      */
 191     bool currentValuesSelected() const
 192     {
 193         return selections_.empty() || currentPosition(0).selected();
 194     }
 195
 196     //! Returns the currently active selection.
 197     const Selection& currentSelection() const
 198     {
 199         GMX_ASSERT(currentSelection_ < ssize(selections_), "Accessing an invalid selection");
 200         return selections_[currentSelection_];
 201     }
 202     //! Returns the `i`th position for the current value.
 203     SelectionPosition currentPosition(int i) const
 204     {
 205         return currentSelection().position(nextPosition_ + i);
 206     }
 207     /*! \brief
 208      * Extracts all coordinates corresponding to the current value.
 209      *
 210      * \param[out] x  Array to which the positions are extracted.
 211      *
 212      * \p x should contain at minimum the number of positions per value
 213      * passed to the constructor.
 214      */
 215     void getCurrentPositions(rvec x[]) const
 216     {
 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)
 221         {
 222             copy_rvec(currentPosition(i).x(), x[i]);
 223         }
 224     }
 225
 226 private:
 227     const SelectionList& selections_;
 228     const int            posCountPerValue_;
 229     int                  currentSelection_;
 230     int                  nextPosition_;
 231
 232     GMX_DISALLOW_COPY_AND_ASSIGN(AnglePositionIterator);
 233 };
 234
 235 /********************************************************************
 236  * Actual analysis module
 237  */
 238
 239 //! How to interpret the selections in -group1.
 240 enum class Group1Type : int
 241 {
 242     Angle,
 243     Dihedral,
 244     Vector,
 245     Plane,
 246     Count
 247 };
 248 //! How to interpret the selections in -group2.
 249 enum class Group2Type : int
 250 {
 251     None,
 252     Vector,
 253     Plane,
 254     TimeZero,
 255     Z,
 256     SphereNormal,
 257     Count
 258 };
 259 //! String values corresponding to Group1Type.
 260 const EnumerationArray<Group1Type, const char*> c_group1TypeEnumNames = {
 261     { "angle", "dihedral", "vector", "plane" }
 262 };
 263 //! String values corresponding to Group2Type.
 264 const EnumerationArray<Group2Type, const char*> c_group2TypeEnumNames = {
 265     { "none", "vector", "plane", "t0", "z", "sphnorm" }
 266 };
 267
 268 class Angle : public TrajectoryAnalysisModule
 269 {
 270 public:
 271     Angle();
 272
 273     void initOptions(IOptionsContainer* options, TrajectoryAnalysisSettings* settings) override;
 274     void optionsFinished(TrajectoryAnalysisSettings* settings) override;
 275     void initAnalysis(const TrajectoryAnalysisSettings& settings, const TopologyInformation& top) override;
 276
 277     void analyzeFrame(int frnr, const t_trxframe& fr, t_pbc* pbc, TrajectoryAnalysisModuleData* pdata) override;
 278
 279     void finishAnalysis(int nframes) override;
 280     void writeOutput() override;
 281
 282 private:
 283     void initFromSelections(const SelectionList& sel1, const SelectionList& sel2);
 284     void checkSelections(const SelectionList& sel1, const SelectionList& sel2) const;
 285
 286     SelectionList        sel1_;
 287     SelectionList        sel2_;
 288     SelectionOptionInfo* sel1info_;
 289     SelectionOptionInfo* sel2info_;
 290     std::string          fnAverage_;
 291     std::string          fnAll_;
 292     std::string          fnHistogram_;
 293
 294     Group1Type g1type_;
 295     Group2Type g2type_;
 296     double     binWidth_;
 297
 298     AnalysisData                             angles_;
 299     AnalysisDataFrameAverageModulePointer    averageModule_;
 300     AnalysisDataSimpleHistogramModulePointer histogramModule_;
 301
 302     std::vector<int>               angleCount_;
 303     int                            natoms1_;
 304     int                            natoms2_;
 305     std::vector<std::vector<RVec>> vt0_;
 306
 307     // Copy and assign disallowed by base.
 308 };
 309
 310 Angle::Angle() :
 311     sel1info_(nullptr),
 312     sel2info_(nullptr),
 313     g1type_(Group1Type::Angle),
 314     g2type_(Group2Type::None),
 315     binWidth_(1.0),
 316     natoms1_(0),
 317     natoms2_(0)
 318 {
 319     averageModule_ = std::make_unique<AnalysisDataFrameAverageModule>();
 320     angles_.addModule(averageModule_);
 321     histogramModule_ = std::make_unique<AnalysisDataSimpleHistogramModule>();
 322     angles_.addModule(histogramModule_);
 323
 324     registerAnalysisDataset(&angles_, "angle");
 325     registerBasicDataset(averageModule_.get(), "average");
 326     registerBasicDataset(&histogramModule_->averager(), "histogram");
 327 }
 328
 329
 330 void Angle::initOptions(IOptionsContainer* options, TrajectoryAnalysisSettings* settings)
 331 {
 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",
 340         "positions.[PAR]",
 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",
 362         "selection.[PAR]",
 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",
 373         "for each frame.",
 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]."
 379     };
 380
 381     settings->setHelpText(desc);
 382
 383     options->addOption(FileNameOption("oav")
 384                                .filetype(OptionFileType::Plot)
 385                                .outputFile()
 386                                .store(&fnAverage_)
 387                                .defaultBasename("angaver")
 388                                .description("Average angles as a function of time"));
 389     options->addOption(FileNameOption("oall")
 390                                .filetype(OptionFileType::Plot)
 391                                .outputFile()
 392                                .store(&fnAll_)
 393                                .defaultBasename("angles")
 394                                .description("All angles as a function of time"));
 395     options->addOption(FileNameOption("oh")
 396                                .filetype(OptionFileType::Plot)
 397                                .outputFile()
 398                                .store(&fnHistogram_)
 399                                .defaultBasename("anghist")
 400                                .description("Histogram of the angles"));
 401
 402     options->addOption(EnumOption<Group1Type>("g1")
 403                                .enumValue(c_group1TypeEnumNames)
 404                                .store(&g1type_)
 405                                .description("Type of analysis/first vector group"));
 406     options->addOption(
 407             EnumOption<Group2Type>("g2").enumValue(c_group2TypeEnumNames).store(&g2type_).description("Type of second vector group"));
 408     options->addOption(
 409             DoubleOption("binw").store(&binWidth_).description("Binwidth for -oh in degrees"));
 410
 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"));
 417 }
 418
 419
 420 void Angle::optionsFinished(TrajectoryAnalysisSettings* /* settings */)
 421 {
 422     const bool bSingle = (g1type_ == Group1Type::Angle || g1type_ == Group1Type::Dihedral);
 423
 424     if (bSingle && g2type_ != Group2Type::None)
 425     {
 426         GMX_THROW(
 427                 InconsistentInputError("Cannot use a second group (-g2) with "
 428                                        "-g1 angle or dihedral"));
 429     }
 430     if (bSingle && sel2info_->isSet())
 431     {
 432         GMX_THROW(
 433                 InconsistentInputError("Cannot provide a second selection "
 434                                        "(-group2) with -g1 angle or dihedral"));
 435     }
 436     if (!bSingle && g2type_ == Group2Type::None)
 437     {
 438         GMX_THROW(
 439                 InconsistentInputError("Should specify a second group (-g2) "
 440                                        "if the first group is not an angle or a dihedral"));
 441     }
 442
 443     // Set up the number of positions per angle.
 444     switch (g1type_)
 445     {
 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"));
 451     }
 452     switch (g2type_)
 453     {
 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"));
 461     }
 462     if (natoms2_ == 0 && sel2info_->isSet())
 463     {
 464         GMX_THROW(InconsistentInputError(
 465                 "Cannot provide a second selection (-group2) with -g2 t0 or z"));
 466     }
 467     // TODO: If bSingle is not set, the second selection option should be
 468     // required.
 469 }
 470
 471
 472 void Angle::initFromSelections(const SelectionList& sel1, const SelectionList& sel2)
 473 {
 474     const int  angleGroups         = std::max(sel1.size(), sel2.size());
 475     const bool bHasSecondSelection = natoms2_ > 0;
 476
 477     if (bHasSecondSelection && sel1.size() != sel2.size() && std::min(sel1.size(), sel2.size()) != 1)
 478     {
 479         GMX_THROW(InconsistentInputError(
 480                 "-group1 and -group2 should specify the same number of selections"));
 481     }
 482
 483     AnglePositionIterator iter1(sel1, natoms1_);
 484     AnglePositionIterator iter2(sel2, natoms2_);
 485     for (int g = 0; g < angleGroups; ++g, iter1.nextGroup(), iter2.nextGroup())
 486     {
 487         const int posCount1 = iter1.currentSelection().posCount();
 488         if (natoms1_ > 1 && posCount1 % natoms1_ != 0)
 489         {
 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),
 493                     natoms1_)));
 494         }
 495         const int angleCount1 = posCount1 / natoms1_;
 496         int       angleCount  = angleCount1;
 497
 498         if (bHasSecondSelection)
 499         {
 500             const int posCount2 = iter2.currentSelection().posCount();
 501             if (natoms2_ > 1 && posCount2 % natoms2_ != 0)
 502             {
 503                 GMX_THROW(InconsistentInputError(
 504                         formatString("Number of positions in selection %d in the second group not "
 505                                      "divisible by %d",
 506                                      static_cast<int>(g + 1),
 507                                      natoms2_)));
 508             }
 509             if (g2type_ == Group2Type::SphereNormal && posCount2 != 1)
 510             {
 511                 GMX_THROW(InconsistentInputError(
 512                         "The second group should contain a single position with -g2 sphnorm"));
 513             }
 514
 515             const int angleCount2 = posCount2 / natoms2_;
 516             angleCount            = std::max(angleCount1, angleCount2);
 517             if (angleCount1 != angleCount2 && std::min(angleCount1, angleCount2) != 1)
 518             {
 519                 GMX_THROW(InconsistentInputError(
 520                         "Number of vectors defined by the two groups are not the same"));
 521             }
 522         }
 523         angleCount_.push_back(angleCount);
 524     }
 525 }
 526
 527
 528 void Angle::checkSelections(const SelectionList& sel1, const SelectionList& sel2) const
 529 {
 530     AnglePositionIterator iter1(sel1, natoms1_);
 531     AnglePositionIterator iter2(sel2, natoms2_);
 532     for (size_t g = 0; g < angleCount_.size(); ++g, iter1.nextGroup(), iter2.nextGroup())
 533     {
 534         if (iter1.isDynamic() || (iter2.hasValue() && iter2.isDynamic()))
 535         {
 536             for (int n = 0; n < angleCount_[g]; ++n, iter1.nextValue(), iter2.nextValue())
 537             {
 538                 bool bOk = true;
 539                 if (!iter1.allValuesConsistentlySelected())
 540                 {
 541                     bOk = false;
 542                 }
 543                 if (!iter2.allValuesConsistentlySelected())
 544                 {
 545                     bOk = false;
 546                 }
 547                 if (angleCount_[g] > 1)
 548                 {
 549                     if (iter1.hasSingleValue() && !iter1.currentValuesSelected())
 550                     {
 551                         bOk = false;
 552                     }
 553                     if (iter2.hasSingleValue() && !iter2.currentValuesSelected())
 554                     {
 555                         bOk = false;
 556                     }
 557                 }
 558                 if (iter2.hasValue()
 559                     && (angleCount_[g] == 1 || (!iter1.hasSingleValue() && !iter2.hasSingleValue()))
 560                     && iter1.currentValuesSelected() != iter2.currentValuesSelected())
 561                 {
 562                     bOk = false;
 563                 }
 564                 if (!bOk)
 565                 {
 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));
 571                 }
 572             }
 573         }
 574     }
 575 }
 576
 577
 578 void Angle::initAnalysis(const TrajectoryAnalysisSettings& settings, const TopologyInformation& /* top */)
 579 {
 580     initFromSelections(sel1_, sel2_);
 581
 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)
 585     {
 586         angles_.setColumnCount(i, angleCount_[i]);
 587     }
 588     double histogramMin = (g1type_ == Group1Type::Dihedral ? -180.0 : 0);
 589     histogramModule_->init(histogramFromRange(histogramMin, 180.0).binWidth(binWidth_).includeAll());
 590
 591     if (g2type_ == Group2Type::TimeZero)
 592     {
 593         vt0_.resize(sel1_.size());
 594         for (size_t g = 0; g < sel1_.size(); ++g)
 595         {
 596             vt0_[g].resize(sel1_[g].posCount() / natoms1_);
 597         }
 598     }
 599
 600     if (!fnAverage_.empty())
 601     {
 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)
 610         {
 611             plotm->appendLegend(sel1_[g].name());
 612         }
 613         averageModule_->addModule(plotm);
 614     }
 615
 616     if (!fnAll_.empty())
 617     {
 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);
 625     }
 626
 627     if (!fnHistogram_.empty())
 628     {
 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)
 637         {
 638             plotm->appendLegend(sel1_[g].name());
 639         }
 640         histogramModule_->averager().addModule(plotm);
 641     }
 642 }
 643
 644
 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)
 647 {
 648     switch (natoms)
 649     {
 650         case 2:
 651             if (pbc)
 652             {
 653                 pbc_dx(pbc, x[1], x[0], xout);
 654             }
 655             else
 656             {
 657                 rvec_sub(x[1], x[0], xout);
 658             }
 659             svmul(0.5, xout, cout);
 660             rvec_add(x[0], cout, cout);
 661             break;
 662         case 3:
 663         {
 664             rvec v1, v2;
 665             if (pbc)
 666             {
 667                 pbc_dx(pbc, x[1], x[0], v1);
 668                 pbc_dx(pbc, x[2], x[0], v2);
 669             }
 670             else
 671             {
 672                 rvec_sub(x[1], x[0], v1);
 673                 rvec_sub(x[2], x[0], v2);
 674             }
 675             cprod(v1, v2, xout);
 676             rvec_add(x[0], x[1], cout);
 677             rvec_add(cout, x[2], cout);
 678             svmul(1.0 / 3.0, cout, cout);
 679             break;
 680         }
 681         default: GMX_RELEASE_ASSERT(false, "Incorrectly initialized number of atoms");
 682     }
 683 }
 684
 685
 686 void Angle::analyzeFrame(int frnr, const t_trxframe& fr, t_pbc* pbc, TrajectoryAnalysisModuleData* pdata)
 687 {
 688     AnalysisDataHandle   dh   = pdata->dataHandle(angles_);
 689     const SelectionList& sel1 = pdata->parallelSelections(sel1_);
 690     const SelectionList& sel2 = pdata->parallelSelections(sel2_);
 691
 692     checkSelections(sel1, sel2);
 693
 694     dh.startFrame(frnr, fr.time);
 695
 696     AnglePositionIterator iter1(sel1, natoms1_);
 697     AnglePositionIterator iter2(sel2, natoms2_);
 698     for (size_t g = 0; g < angleCount_.size(); ++g, iter1.nextGroup(), iter2.nextGroup())
 699     {
 700         rvec v1, v2;
 701         rvec c1, c2;
 702
 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.
 708
 709         clear_rvec(v2);
 710         clear_rvec(c2);
 711
 712         switch (g2type_)
 713         {
 714             case Group2Type::Z: v2[ZZ] = 1.0; break;
 715             case Group2Type::SphereNormal: copy_rvec(sel2_[g].position(0).x(), c2); break;
 716             default:
 717                 // do nothing
 718                 break;
 719         }
 720
 721         dh.selectDataSet(g);
 722         for (int n = 0; n < angleCount_[g]; ++n, iter1.nextValue(), iter2.nextValue())
 723         {
 724             rvec x[4];
 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.
 728             clear_rvecs(4, x);
 729
 730             real angle = 0;
 731             // checkSelections() ensures that this reflects all the involved
 732             // positions.
 733             const bool bPresent = iter1.currentValuesSelected() && iter2.currentValuesSelected();
 734             iter1.getCurrentPositions(x);
 735             switch (g1type_)
 736             {
 737                 case Group1Type::Angle:
 738                     if (pbc)
 739                     {
 740                         pbc_dx(pbc, x[0], x[1], v1);
 741                         pbc_dx(pbc, x[2], x[1], v2);
 742                     }
 743                     else
 744                     {
 745                         rvec_sub(x[0], x[1], v1);
 746                         rvec_sub(x[2], x[1], v2);
 747                     }
 748                     angle = gmx_angle(v1, v2);
 749                     break;
 750                 case Group1Type::Dihedral:
 751                 {
 752                     rvec dx[3];
 753                     if (pbc)
 754                     {
 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]);
 758                     }
 759                     else
 760                     {
 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]);
 764                     }
 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);
 769                     if (ipr < 0)
 770                     {
 771                         angle = -angle;
 772                     }
 773                     break;
 774                 }
 775                 case Group1Type::Vector:
 776                 case Group1Type::Plane:
 777                     calc_vec(natoms1_, x, pbc, v1, c1);
 778                     switch (g2type_)
 779                     {
 780                         case Group2Type::Vector:
 781                         case Group2Type::Plane:
 782                             iter2.getCurrentPositions(x);
 783                             calc_vec(natoms2_, x, pbc, v2, c2);
 784                             break;
 785                         case Group2Type::TimeZero:
 786                             // FIXME: This is not parallelizable.
 787                             if (frnr == 0)
 788                             {
 789                                 copy_rvec(v1, vt0_[g][n]);
 790                             }
 791                             copy_rvec(vt0_[g][n], v2);
 792                             break;
 793                         case Group2Type::Z: c1[XX] = c1[YY] = 0.0; break;
 794                         case Group2Type::SphereNormal:
 795                             if (pbc)
 796                             {
 797                                 pbc_dx(pbc, c1, c2, v2);
 798                             }
 799                             else
 800                             {
 801                                 rvec_sub(c1, c2, v2);
 802                             }
 803                             break;
 804                         default: GMX_THROW(InternalError("invalid -g2 value"));
 805                     }
 806                     angle = gmx_angle(v1, v2);
 807                     break;
 808                 default: GMX_THROW(InternalError("invalid -g1 value"));
 809             }
 810             dh.setPoint(n, angle * gmx::c_rad2Deg, bPresent);
 811         }
 812     }
 813     dh.finishFrame();
 814 }
 815
 816
 817 void Angle::finishAnalysis(int /*nframes*/)
 818 {
 819     AbstractAverageHistogram& averageHistogram = histogramModule_->averager();
 820     averageHistogram.normalizeProbability();
 821     averageHistogram.done();
 822 }
 823
 824
 825 void Angle::writeOutput() {}
 826
 827 } // namespace
 828
 829 const char AngleInfo::name[]             = "gangle";
 830 const char AngleInfo::shortDescription[] = "Calculate angles";
 831
 832 TrajectoryAnalysisModulePointer AngleInfo::create()
 833 {
 834     return TrajectoryAnalysisModulePointer(new Angle);
 835 }
 836
 837 } // namespace analysismodules
 838
 839 } // namespace gmx