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