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