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