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/vec.h"
  58 #include "gromacs/options/basicoptions.h"
  59 #include "gromacs/options/filenameoption.h"
  60 #include "gromacs/options/ioptionscontainer.h"
  61 #include "gromacs/pbcutil/pbc.h"
  62 #include "gromacs/selection/selection.h"
  63 #include "gromacs/selection/selectionoption.h"
  64 #include "gromacs/trajectory/trajectoryframe.h"
  65 #include "gromacs/trajectoryanalysis/analysissettings.h"
  66 #include "gromacs/utility/arrayref.h"
  67 #include "gromacs/utility/classhelpers.h"
  68 #include "gromacs/utility/enumerationhelpers.h"
  69 #include "gromacs/utility/exceptions.h"
  70 #include "gromacs/utility/gmxassert.h"
  71 #include "gromacs/utility/stringutil.h"
  72
  73 namespace gmx
  74 {
  75
  76 namespace analysismodules
  77 {
  78
  79 namespace
  80 {
  81
  82 /********************************************************************
  83  * Helper classes
  84  */
  85
  86 /*! \brief
  87  * Helper to encapsulate logic for looping over input selections.
  88  *
  89  * This class provides two-dimensional iteration:
  90  *  - Over _angle groups_, corresponding to an input selection.  If the input
  91  *    selection list contains a single selection, that selection gets used
  92  *    for all angle groups.
  93  *  - Within a group, over _values_, each consisting of a fixed number of
  94  *    selection positions.  If there is only a single value within a selection,
  95  *    that value is returned over and over again.
  96  * This transparently provides the semantics of using a single selection/vector
  97  * to compute angles against multiple selections/vectors as described in the
  98  * tool help text.
  99  *
 100  * This class isn't perferctly self-contained and requires the caller to know
 101  * some of the internals to use it properly, but it serves its purpose for this
 102  * single analysis tool by simplifying the loops.
 103  * Some methods have also been tailored to allow the caller to use it a bit
 104  * more easily.
 105  *
 106  * \ingroup module_trajectoryanalysis
 107  */
 108 class AnglePositionIterator
 109 {
 110 public:
 111     /*! \brief
 112      * Creates an iterator to loop over input selection positions.
 113      *
 114      * \param[in] selections       List of selections.
 115      * \param[in] posCountPerValue Number of selection positions that
 116      *     constitute a single value for the iteration.
 117      *
 118      * If \p selections is empty, and/or \p posCountPerValue is zero, the
 119      * iterator can still be advanced and hasValue()/hasSingleValue()
 120      * called, but values cannot be accessed.
 121      */
 122     AnglePositionIterator(const SelectionList& selections, int posCountPerValue) :
 123         selections_(selections),
 124         posCountPerValue_(posCountPerValue),
 125         currentSelection_(0),
 126         nextPosition_(0)
 127     {
 128     }
 129
 130     //! Advances the iterator to the next group of angles.
 131     void nextGroup()
 132     {
 133         if (selections_.size() > 1)
 134         {
 135             ++currentSelection_;
 136         }
 137         nextPosition_ = 0;
 138     }
 139     //! Advances the iterator to the next angle in the current group.
 140     void nextValue()
 141     {
 142         if (!hasSingleValue())
 143         {
 144             nextPosition_ += posCountPerValue_;
 145         }
 146     }
 147
 148     /*! \brief
 149      * Returns whether this iterator represents any values.
 150      *
 151      * If the return value is `false`, only nextGroup(), nextValue() and
 152      * hasSingleValue() are allowed to be called.
 153      */
 154     bool hasValue() const { return !selections_.empty(); }
 155     /*! \brief
 156      * Returns whether the current selection only contains a single value.
 157      *
 158      * Returns `false` if hasValue() returns false, which allows cutting
 159      * some corners in consistency checks.
 160      */
 161     bool hasSingleValue() const
 162     {
 163         return hasValue() && currentSelection().posCount() == posCountPerValue_;
 164     }
 165     //! Returns whether the current selection is dynamic.
 166     bool isDynamic() const { return currentSelection().isDynamic(); }
 167     /*! \brief
 168      * Returns whether positions in the current value are either all
 169      * selected or all unselected.
 170      */
 171     bool allValuesConsistentlySelected() const
 172     {
 173         if (posCountPerValue_ <= 1)
 174         {
 175             return true;
 176         }
 177         const bool bSelected = currentPosition(0).selected();
 178         for (int i = 1; i < posCountPerValue_; ++i)
 179         {
 180             if (currentPosition(i).selected() != bSelected)
 181             {
 182                 return false;
 183             }
 184         }
 185         return true;
 186     }
 187     /*! \brief
 188      * Returns whether positions in the current value are selected.
 189      *
 190      * Only works reliably if allValuesConsistentlySelected() returns
 191      * `true`.
 192      */
 193     bool currentValuesSelected() const
 194     {
 195         return selections_.empty() || currentPosition(0).selected();
 196     }
 197
 198     //! Returns the currently active selection.
 199     const Selection& currentSelection() const
 200     {
 201         GMX_ASSERT(currentSelection_ < ssize(selections_), "Accessing an invalid selection");
 202         return selections_[currentSelection_];
 203     }
 204     //! Returns the `i`th position for the current value.
 205     SelectionPosition currentPosition(int i) const
 206     {
 207         return currentSelection().position(nextPosition_ + i);
 208     }
 209     /*! \brief
 210      * Extracts all coordinates corresponding to the current value.
 211      *
 212      * \param[out] x  Array to which the positions are extracted.
 213      *
 214      * \p x should contain at minimum the number of positions per value
 215      * passed to the constructor.
 216      */
 217     void getCurrentPositions(rvec x[]) const
 218     {
 219         GMX_ASSERT(posCountPerValue_ > 0, "Accessing positions for an invalid angle type");
 220         GMX_ASSERT(nextPosition_ + posCountPerValue_ <= currentSelection().posCount(),
 221                    "Accessing an invalid position");
 222         for (int i = 0; i < posCountPerValue_; ++i)
 223         {
 224             copy_rvec(currentPosition(i).x(), x[i]);
 225         }
 226     }
 227
 228 private:
 229     const SelectionList& selections_;
 230     const int            posCountPerValue_;
 231     int                  currentSelection_;
 232     int                  nextPosition_;
 233
 234     GMX_DISALLOW_COPY_AND_ASSIGN(AnglePositionIterator);
 235 };
 236
 237 /********************************************************************
 238  * Actual analysis module
 239  */
 240
 241 //! How to interpret the selections in -group1.
 242 enum class Group1Type : int
 243 {
 244     Angle,
 245     Dihedral,
 246     Vector,
 247     Plane,
 248     Count
 249 };
 250 //! How to interpret the selections in -group2.
 251 enum class Group2Type : int
 252 {
 253     None,
 254     Vector,
 255     Plane,
 256     TimeZero,
 257     Z,
 258     SphereNormal,
 259     Count
 260 };
 261 //! String values corresponding to Group1Type.
 262 const EnumerationArray<Group1Type, const char*> c_group1TypeEnumNames = {
 263     { "angle", "dihedral", "vector", "plane" }
 264 };
 265 //! String values corresponding to Group2Type.
 266 const EnumerationArray<Group2Type, const char*> c_group2TypeEnumNames = {
 267     { "none", "vector", "plane", "t0", "z", "sphnorm" }
 268 };
 269
 270 class Angle : public TrajectoryAnalysisModule
 271 {
 272 public:
 273     Angle();
 274
 275     void initOptions(IOptionsContainer* options, TrajectoryAnalysisSettings* settings) override;
 276     void optionsFinished(TrajectoryAnalysisSettings* settings) override;
 277     void initAnalysis(const TrajectoryAnalysisSettings& settings, const TopologyInformation& top) override;
 278
 279     void analyzeFrame(int frnr, const t_trxframe& fr, t_pbc* pbc, TrajectoryAnalysisModuleData* pdata) override;
 280
 281     void finishAnalysis(int nframes) override;
 282     void writeOutput() override;
 283
 284 private:
 285     void initFromSelections(const SelectionList& sel1, const SelectionList& sel2);
 286     void checkSelections(const SelectionList& sel1, const SelectionList& sel2) const;
 287
 288     SelectionList        sel1_;
 289     SelectionList        sel2_;
 290     SelectionOptionInfo* sel1info_;
 291     SelectionOptionInfo* sel2info_;
 292     std::string          fnAverage_;
 293     std::string          fnAll_;
 294     std::string          fnHistogram_;
 295
 296     Group1Type g1type_;
 297     Group2Type g2type_;
 298     double     binWidth_;
 299
 300     AnalysisData                             angles_;
 301     AnalysisDataFrameAverageModulePointer    averageModule_;
 302     AnalysisDataSimpleHistogramModulePointer histogramModule_;
 303
 304     std::vector<int>               angleCount_;
 305     int                            natoms1_;
 306     int                            natoms2_;
 307     std::vector<std::vector<RVec>> vt0_;
 308
 309     // Copy and assign disallowed by base.
 310 };
 311
 312 Angle::Angle() :
 313     sel1info_(nullptr),
 314     sel2info_(nullptr),
 315     g1type_(Group1Type::Angle),
 316     g2type_(Group2Type::None),
 317     binWidth_(1.0),
 318     natoms1_(0),
 319     natoms2_(0)
 320 {
 321     averageModule_ = std::make_unique<AnalysisDataFrameAverageModule>();
 322     angles_.addModule(averageModule_);
 323     histogramModule_ = std::make_unique<AnalysisDataSimpleHistogramModule>();
 324     angles_.addModule(histogramModule_);
 325
 326     registerAnalysisDataset(&angles_, "angle");
 327     registerBasicDataset(averageModule_.get(), "average");
 328     registerBasicDataset(&histogramModule_->averager(), "histogram");
 329 }
 330
 331
 332 void Angle::initOptions(IOptionsContainer* options, TrajectoryAnalysisSettings* settings)
 333 {
 334     static const char* const desc[] = {
 335         "[THISMODULE] computes different types of angles between vectors.",
 336         "It supports both vectors defined by two positions and normals of",
 337         "planes defined by three positions.",
 338         "The z axis or the local normal of a sphere can also be used as",
 339         "one of the vectors.",
 340         "There are also convenience options 'angle' and 'dihedral' for",
 341         "calculating bond angles and dihedrals defined by three/four",
 342         "positions.[PAR]",
 343         "The type of the angle is specified with [TT]-g1[tt] and [TT]-g2[tt].",
 344         "If [TT]-g1[tt] is [TT]angle[tt] or [TT]dihedral[tt], [TT]-g2[tt]",
 345         "should not be specified.",
 346         "In this case, [TT]-group1[tt] should specify one or more selections,",
 347         "and each should contain triplets or quartets of positions that define",
 348         "the angles to be calculated.[PAR]",
 349         "If [TT]-g1[tt] is [TT]vector[tt] or [TT]plane[tt], [TT]-group1[tt]",
 350         "should specify selections that contain either pairs ([TT]vector[tt])",
 351         "or triplets ([TT]plane[tt]) of positions. For vectors, the positions",
 352         "set the endpoints of the vector, and for planes, the three positions",
 353         "are used to calculate the normal of the plane. In both cases,",
 354         "[TT]-g2[tt] specifies the other vector to use (see below).[PAR]",
 355         "With [TT]-g2 vector[tt] or [TT]-g2 plane[tt], [TT]-group2[tt] should",
 356         "specify another set of vectors. [TT]-group1[tt] and [TT]-group2[tt]",
 357         "should specify the same number of selections. It is also allowed to",
 358         "only have a single selection for one of the options, in which case",
 359         "the same selection is used with each selection in the other group.",
 360         "Similarly, for each selection in [TT]-group1[tt], the corresponding",
 361         "selection in [TT]-group2[tt] should specify the same number of",
 362         "vectors or a single vector. In the latter case, the angle is",
 363         "calculated between that single vector and each vector from the other",
 364         "selection.[PAR]",
 365         "With [TT]-g2 sphnorm[tt], each selection in [TT]-group2[tt] should",
 366         "specify a single position that is the center of the sphere.",
 367         "The second vector is calculated as the vector from the center to the",
 368         "midpoint of the positions specified by [TT]-group1[tt].[PAR]",
 369         "With [TT]-g2 z[tt], [TT]-group2[tt] is not necessary, and angles",
 370         "between the first vectors and the positive Z axis are calculated.[PAR]",
 371         "With [TT]-g2 t0[tt], [TT]-group2[tt] is not necessary, and angles",
 372         "are calculated from the vectors as they are in the first frame.[PAR]",
 373         "There are three options for output:",
 374         "[TT]-oav[tt] writes an xvg file with the time and the average angle",
 375         "for each frame.",
 376         "[TT]-oall[tt] writes all the individual angles.",
 377         "[TT]-oh[tt] writes a histogram of the angles. The bin width can be",
 378         "set with [TT]-binw[tt].",
 379         "For [TT]-oav[tt] and [TT]-oh[tt], separate average/histogram is",
 380         "computed for each selection in [TT]-group1[tt]."
 381     };
 382
 383     settings->setHelpText(desc);
 384
 385     options->addOption(FileNameOption("oav")
 386                                .filetype(eftPlot)
 387                                .outputFile()
 388                                .store(&fnAverage_)
 389                                .defaultBasename("angaver")
 390                                .description("Average angles as a function of time"));
 391     options->addOption(FileNameOption("oall")
 392                                .filetype(eftPlot)
 393                                .outputFile()
 394                                .store(&fnAll_)
 395                                .defaultBasename("angles")
 396                                .description("All angles as a function of time"));
 397     options->addOption(FileNameOption("oh")
 398                                .filetype(eftPlot)
 399                                .outputFile()
 400                                .store(&fnHistogram_)
 401                                .defaultBasename("anghist")
 402                                .description("Histogram of the angles"));
 403
 404     options->addOption(EnumOption<Group1Type>("g1")
 405                                .enumValue(c_group1TypeEnumNames)
 406                                .store(&g1type_)
 407                                .description("Type of analysis/first vector group"));
 408     options->addOption(
 409             EnumOption<Group2Type>("g2").enumValue(c_group2TypeEnumNames).store(&g2type_).description("Type of second vector group"));
 410     options->addOption(
 411             DoubleOption("binw").store(&binWidth_).description("Binwidth for -oh in degrees"));
 412
 413     sel1info_ = options->addOption(
 414             SelectionOption("group1").required().dynamicMask().storeVector(&sel1_).multiValue().description(
 415                     "First analysis/vector selection"));
 416     sel2info_ = options->addOption(
 417             SelectionOption("group2").dynamicMask().storeVector(&sel2_).multiValue().description(
 418                     "Second analysis/vector selection"));
 419 }
 420
 421
 422 void Angle::optionsFinished(TrajectoryAnalysisSettings* /* settings */)
 423 {
 424     const bool bSingle = (g1type_ == Group1Type::Angle || g1type_ == Group1Type::Dihedral);
 425
 426     if (bSingle && g2type_ != Group2Type::None)
 427     {
 428         GMX_THROW(
 429                 InconsistentInputError("Cannot use a second group (-g2) with "
 430                                        "-g1 angle or dihedral"));
 431     }
 432     if (bSingle && sel2info_->isSet())
 433     {
 434         GMX_THROW(
 435                 InconsistentInputError("Cannot provide a second selection "
 436                                        "(-group2) with -g1 angle or dihedral"));
 437     }
 438     if (!bSingle && g2type_ == Group2Type::None)
 439     {
 440         GMX_THROW(
 441                 InconsistentInputError("Should specify a second group (-g2) "
 442                                        "if the first group is not an angle or a dihedral"));
 443     }
 444
 445     // Set up the number of positions per angle.
 446     switch (g1type_)
 447     {
 448         case Group1Type::Angle: natoms1_ = 3; break;
 449         case Group1Type::Dihedral: natoms1_ = 4; break;
 450         case Group1Type::Vector: natoms1_ = 2; break;
 451         case Group1Type::Plane: natoms1_ = 3; break;
 452         default: GMX_THROW(InternalError("invalid -g1 value"));
 453     }
 454     switch (g2type_)
 455     {
 456         case Group2Type::None: natoms2_ = 0; break;
 457         case Group2Type::Vector: natoms2_ = 2; break;
 458         case Group2Type::Plane: natoms2_ = 3; break;
 459         case Group2Type::TimeZero: // Intended to fall through
 460         case Group2Type::Z: natoms2_ = 0; break;
 461         case Group2Type::SphereNormal: natoms2_ = 1; break;
 462         default: GMX_THROW(InternalError("invalid -g2 value"));
 463     }
 464     if (natoms2_ == 0 && sel2info_->isSet())
 465     {
 466         GMX_THROW(InconsistentInputError(
 467                 "Cannot provide a second selection (-group2) with -g2 t0 or z"));
 468     }
 469     // TODO: If bSingle is not set, the second selection option should be
 470     // required.
 471 }
 472
 473
 474 void Angle::initFromSelections(const SelectionList& sel1, const SelectionList& sel2)
 475 {
 476     const int  angleGroups         = std::max(sel1.size(), sel2.size());
 477     const bool bHasSecondSelection = natoms2_ > 0;
 478
 479     if (bHasSecondSelection && sel1.size() != sel2.size() && std::min(sel1.size(), sel2.size()) != 1)
 480     {
 481         GMX_THROW(InconsistentInputError(
 482                 "-group1 and -group2 should specify the same number of selections"));
 483     }
 484
 485     AnglePositionIterator iter1(sel1, natoms1_);
 486     AnglePositionIterator iter2(sel2, natoms2_);
 487     for (int g = 0; g < angleGroups; ++g, iter1.nextGroup(), iter2.nextGroup())
 488     {
 489         const int posCount1 = iter1.currentSelection().posCount();
 490         if (natoms1_ > 1 && posCount1 % natoms1_ != 0)
 491         {
 492             GMX_THROW(InconsistentInputError(formatString(
 493                     "Number of positions in selection %d in the first group not divisible by %d",
 494                     static_cast<int>(g + 1),
 495                     natoms1_)));
 496         }
 497         const int angleCount1 = posCount1 / natoms1_;
 498         int       angleCount  = angleCount1;
 499
 500         if (bHasSecondSelection)
 501         {
 502             const int posCount2 = iter2.currentSelection().posCount();
 503             if (natoms2_ > 1 && posCount2 % natoms2_ != 0)
 504             {
 505                 GMX_THROW(InconsistentInputError(
 506                         formatString("Number of positions in selection %d in the second group not "
 507                                      "divisible by %d",
 508                                      static_cast<int>(g + 1),
 509                                      natoms2_)));
 510             }
 511             if (g2type_ == Group2Type::SphereNormal && posCount2 != 1)
 512             {
 513                 GMX_THROW(InconsistentInputError(
 514                         "The second group should contain a single position with -g2 sphnorm"));
 515             }
 516
 517             const int angleCount2 = posCount2 / natoms2_;
 518             angleCount            = std::max(angleCount1, angleCount2);
 519             if (angleCount1 != angleCount2 && std::min(angleCount1, angleCount2) != 1)
 520             {
 521                 GMX_THROW(InconsistentInputError(
 522                         "Number of vectors defined by the two groups are not the same"));
 523             }
 524         }
 525         angleCount_.push_back(angleCount);
 526     }
 527 }
 528
 529
 530 void Angle::checkSelections(const SelectionList& sel1, const SelectionList& sel2) const
 531 {
 532     AnglePositionIterator iter1(sel1, natoms1_);
 533     AnglePositionIterator iter2(sel2, natoms2_);
 534     for (size_t g = 0; g < angleCount_.size(); ++g, iter1.nextGroup(), iter2.nextGroup())
 535     {
 536         if (iter1.isDynamic() || (iter2.hasValue() && iter2.isDynamic()))
 537         {
 538             for (int n = 0; n < angleCount_[g]; ++n, iter1.nextValue(), iter2.nextValue())
 539             {
 540                 bool bOk = true;
 541                 if (!iter1.allValuesConsistentlySelected())
 542                 {
 543                     bOk = false;
 544                 }
 545                 if (!iter2.allValuesConsistentlySelected())
 546                 {
 547                     bOk = false;
 548                 }
 549                 if (angleCount_[g] > 1)
 550                 {
 551                     if (iter1.hasSingleValue() && !iter1.currentValuesSelected())
 552                     {
 553                         bOk = false;
 554                     }
 555                     if (iter2.hasSingleValue() && !iter2.currentValuesSelected())
 556                     {
 557                         bOk = false;
 558                     }
 559                 }
 560                 if (iter2.hasValue()
 561                     && (angleCount_[g] == 1 || (!iter1.hasSingleValue() && !iter2.hasSingleValue()))
 562                     && iter1.currentValuesSelected() != iter2.currentValuesSelected())
 563                 {
 564                     bOk = false;
 565                 }
 566                 if (!bOk)
 567                 {
 568                     std::string message = formatString(
 569                             "Dynamic selection %d does not select "
 570                             "a consistent set of angles over the frames",
 571                             static_cast<int>(g + 1));
 572                     GMX_THROW(InconsistentInputError(message));
 573                 }
 574             }
 575         }
 576     }
 577 }
 578
 579
 580 void Angle::initAnalysis(const TrajectoryAnalysisSettings& settings, const TopologyInformation& /* top */)
 581 {
 582     initFromSelections(sel1_, sel2_);
 583
 584     // checkSelections() ensures that both selection lists have the same size.
 585     angles_.setDataSetCount(angleCount_.size());
 586     for (size_t i = 0; i < angleCount_.size(); ++i)
 587     {
 588         angles_.setColumnCount(i, angleCount_[i]);
 589     }
 590     double histogramMin = (g1type_ == Group1Type::Dihedral ? -180.0 : 0);
 591     histogramModule_->init(histogramFromRange(histogramMin, 180.0).binWidth(binWidth_).includeAll());
 592
 593     if (g2type_ == Group2Type::TimeZero)
 594     {
 595         vt0_.resize(sel1_.size());
 596         for (size_t g = 0; g < sel1_.size(); ++g)
 597         {
 598             vt0_[g].resize(sel1_[g].posCount() / natoms1_);
 599         }
 600     }
 601
 602     if (!fnAverage_.empty())
 603     {
 604         AnalysisDataPlotModulePointer plotm(new AnalysisDataPlotModule(settings.plotSettings()));
 605         plotm->setFileName(fnAverage_);
 606         plotm->setTitle("Average angle");
 607         plotm->setXAxisIsTime();
 608         plotm->setYLabel("Angle (degrees)");
 609         // TODO: Consider adding information about the second selection,
 610         // and/or a subtitle describing what kind of angle this is.
 611         for (size_t g = 0; g < sel1_.size(); ++g)
 612         {
 613             plotm->appendLegend(sel1_[g].name());
 614         }
 615         averageModule_->addModule(plotm);
 616     }
 617
 618     if (!fnAll_.empty())
 619     {
 620         AnalysisDataPlotModulePointer plotm(new AnalysisDataPlotModule(settings.plotSettings()));
 621         plotm->setFileName(fnAll_);
 622         plotm->setTitle("Angle");
 623         plotm->setXAxisIsTime();
 624         plotm->setYLabel("Angle (degrees)");
 625         // TODO: Add legends? (there can be a massive amount of columns)
 626         angles_.addModule(plotm);
 627     }
 628
 629     if (!fnHistogram_.empty())
 630     {
 631         AnalysisDataPlotModulePointer plotm(new AnalysisDataPlotModule(settings.plotSettings()));
 632         plotm->setFileName(fnHistogram_);
 633         plotm->setTitle("Angle histogram");
 634         plotm->setXLabel("Angle (degrees)");
 635         plotm->setYLabel("Probability");
 636         // TODO: Consider adding information about the second selection,
 637         // and/or a subtitle describing what kind of angle this is.
 638         for (size_t g = 0; g < sel1_.size(); ++g)
 639         {
 640             plotm->appendLegend(sel1_[g].name());
 641         }
 642         histogramModule_->averager().addModule(plotm);
 643     }
 644 }
 645
 646
 647 //! Helper method to calculate a vector from two or three positions.
 648 void calc_vec(int natoms, rvec x[], t_pbc* pbc, rvec xout, rvec cout)
 649 {
 650     switch (natoms)
 651     {
 652         case 2:
 653             if (pbc)
 654             {
 655                 pbc_dx(pbc, x[1], x[0], xout);
 656             }
 657             else
 658             {
 659                 rvec_sub(x[1], x[0], xout);
 660             }
 661             svmul(0.5, xout, cout);
 662             rvec_add(x[0], cout, cout);
 663             break;
 664         case 3:
 665         {
 666             rvec v1, v2;
 667             if (pbc)
 668             {
 669                 pbc_dx(pbc, x[1], x[0], v1);
 670                 pbc_dx(pbc, x[2], x[0], v2);
 671             }
 672             else
 673             {
 674                 rvec_sub(x[1], x[0], v1);
 675                 rvec_sub(x[2], x[0], v2);
 676             }
 677             cprod(v1, v2, xout);
 678             rvec_add(x[0], x[1], cout);
 679             rvec_add(cout, x[2], cout);
 680             svmul(1.0 / 3.0, cout, cout);
 681             break;
 682         }
 683         default: GMX_RELEASE_ASSERT(false, "Incorrectly initialized number of atoms");
 684     }
 685 }
 686
 687
 688 void Angle::analyzeFrame(int frnr, const t_trxframe& fr, t_pbc* pbc, TrajectoryAnalysisModuleData* pdata)
 689 {
 690     AnalysisDataHandle   dh   = pdata->dataHandle(angles_);
 691     const SelectionList& sel1 = TrajectoryAnalysisModuleData::parallelSelections(sel1_);
 692     const SelectionList& sel2 = TrajectoryAnalysisModuleData::parallelSelections(sel2_);
 693
 694     checkSelections(sel1, sel2);
 695
 696     dh.startFrame(frnr, fr.time);
 697
 698     AnglePositionIterator iter1(sel1, natoms1_);
 699     AnglePositionIterator iter2(sel2, natoms2_);
 700     for (size_t g = 0; g < angleCount_.size(); ++g, iter1.nextGroup(), iter2.nextGroup())
 701     {
 702         rvec v1, v2;
 703         rvec c1, c2;
 704
 705         // v2 & c2 are conditionally set in the switch statement below, and conditionally
 706         // used in a different switch statement later. Apparently the clang static analyzer
 707         // thinks there are cases where they can be used uninitialzed (which I cannot find),
 708         // but to avoid trouble if we ever change just one of the switch statements it
 709         // makes sense to clear them outside the first switch.
 710
 711         clear_rvec(v2);
 712         clear_rvec(c2);
 713
 714         switch (g2type_)
 715         {
 716             case Group2Type::Z: v2[ZZ] = 1.0; break;
 717             case Group2Type::SphereNormal: copy_rvec(sel2_[g].position(0).x(), c2); break;
 718             default:
 719                 // do nothing
 720                 break;
 721         }
 722
 723         dh.selectDataSet(g);
 724         for (int n = 0; n < angleCount_[g]; ++n, iter1.nextValue(), iter2.nextValue())
 725         {
 726             rvec x[4];
 727             // x[] will be assigned below based on the number of atoms used to initialize iter1,
 728             // which in turn should correspond perfectly to g1type_ (which determines how many we
 729             // read), but unsurprisingly the static analyzer chokes a bit on that.
 730             clear_rvecs(4, x);
 731
 732             real angle = 0;
 733             // checkSelections() ensures that this reflects all the involved
 734             // positions.
 735             const bool bPresent = iter1.currentValuesSelected() && iter2.currentValuesSelected();
 736             iter1.getCurrentPositions(x);
 737             switch (g1type_)
 738             {
 739                 case Group1Type::Angle:
 740                     if (pbc)
 741                     {
 742                         pbc_dx(pbc, x[0], x[1], v1);
 743                         pbc_dx(pbc, x[2], x[1], v2);
 744                     }
 745                     else
 746                     {
 747                         rvec_sub(x[0], x[1], v1);
 748                         rvec_sub(x[2], x[1], v2);
 749                     }
 750                     angle = gmx_angle(v1, v2);
 751                     break;
 752                 case Group1Type::Dihedral:
 753                 {
 754                     rvec dx[3];
 755                     if (pbc)
 756                     {
 757                         pbc_dx(pbc, x[0], x[1], dx[0]);
 758                         pbc_dx(pbc, x[2], x[1], dx[1]);
 759                         pbc_dx(pbc, x[2], x[3], dx[2]);
 760                     }
 761                     else
 762                     {
 763                         rvec_sub(x[0], x[1], dx[0]);
 764                         rvec_sub(x[2], x[1], dx[1]);
 765                         rvec_sub(x[2], x[3], dx[2]);
 766                     }
 767                     cprod(dx[0], dx[1], v1);
 768                     cprod(dx[1], dx[2], v2);
 769                     angle    = gmx_angle(v1, v2);
 770                     real ipr = iprod(dx[0], v2);
 771                     if (ipr < 0)
 772                     {
 773                         angle = -angle;
 774                     }
 775                     break;
 776                 }
 777                 case Group1Type::Vector:
 778                 case Group1Type::Plane:
 779                     calc_vec(natoms1_, x, pbc, v1, c1);
 780                     switch (g2type_)
 781                     {
 782                         case Group2Type::Vector:
 783                         case Group2Type::Plane:
 784                             iter2.getCurrentPositions(x);
 785                             calc_vec(natoms2_, x, pbc, v2, c2);
 786                             break;
 787                         case Group2Type::TimeZero:
 788                             // FIXME: This is not parallelizable.
 789                             if (frnr == 0)
 790                             {
 791                                 copy_rvec(v1, vt0_[g][n]);
 792                             }
 793                             copy_rvec(vt0_[g][n], v2);
 794                             break;
 795                         case Group2Type::Z: c1[XX] = c1[YY] = 0.0; break;
 796                         case Group2Type::SphereNormal:
 797                             if (pbc)
 798                             {
 799                                 pbc_dx(pbc, c1, c2, v2);
 800                             }
 801                             else
 802                             {
 803                                 rvec_sub(c1, c2, v2);
 804                             }
 805                             break;
 806                         default: GMX_THROW(InternalError("invalid -g2 value"));
 807                     }
 808                     angle = gmx_angle(v1, v2);
 809                     break;
 810                 default: GMX_THROW(InternalError("invalid -g1 value"));
 811             }
 812             dh.setPoint(n, angle * RAD2DEG, bPresent);
 813         }
 814     }
 815     dh.finishFrame();
 816 }
 817
 818
 819 void Angle::finishAnalysis(int /*nframes*/)
 820 {
 821     AbstractAverageHistogram& averageHistogram = histogramModule_->averager();
 822     averageHistogram.normalizeProbability();
 823     averageHistogram.done();
 824 }
 825
 826
 827 void Angle::writeOutput() {}
 828
 829 } // namespace
 830
 831 const char AngleInfo::name[]             = "gangle";
 832 const char AngleInfo::shortDescription[] = "Calculate angles";
 833
 834 TrajectoryAnalysisModulePointer AngleInfo::create()
 835 {
 836     return TrajectoryAnalysisModulePointer(new Angle);
 837 }
 838
 839 } // namespace analysismodules
 840
 841 } // namespace gmx