src/gromacs/trajectoryanalysis/modules/angle.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2011,2012,2013,2014,2015,2016,2017, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 /*! \internal \file
  36  * \brief
  37  * Implements gmx::analysismodules::Angle.
  38  *
  39  * \author Teemu Murtola <teemu.murtola@gmail.com>
  40  * \ingroup module_trajectoryanalysis
  41  */
  42 #include "gmxpre.h"
  43
  44 #include "angle.h"
  45
  46 #include <algorithm>
  47 #include <string>
  48 #include <vector>
  49
  50 #include "gromacs/analysisdata/analysisdata.h"
  51 #include "gromacs/analysisdata/modules/average.h"
  52 #include "gromacs/analysisdata/modules/histogram.h"
  53 #include "gromacs/analysisdata/modules/plot.h"
  54 #include "gromacs/math/units.h"
  55 #include "gromacs/math/vec.h"
  56 #include "gromacs/options/basicoptions.h"
  57 #include "gromacs/options/filenameoption.h"
  58 #include "gromacs/options/ioptionscontainer.h"
  59 #include "gromacs/pbcutil/pbc.h"
  60 #include "gromacs/selection/selection.h"
  61 #include "gromacs/selection/selectionoption.h"
  62 #include "gromacs/trajectory/trajectoryframe.h"
  63 #include "gromacs/trajectoryanalysis/analysissettings.h"
  64 #include "gromacs/utility/arrayref.h"
  65 #include "gromacs/utility/exceptions.h"
  66 #include "gromacs/utility/gmxassert.h"
  67 #include "gromacs/utility/stringutil.h"
  68
  69 namespace gmx
  70 {
  71
  72 namespace analysismodules
  73 {
  74
  75 namespace
  76 {
  77
  78 /********************************************************************
  79  * Helper classes
  80  */
  81
  82 /*! \brief
  83  * Helper to encapsulate logic for looping over input selections.
  84  *
  85  * This class provides two-dimensional iteration:
  86  *  - Over _angle groups_, corresponding to an input selection.  If the input
  87  *    selection list contains a single selection, that selection gets used
  88  *    for all angle groups.
  89  *  - Within a group, over _values_, each consisting of a fixed number of
  90  *    selection positions.  If there is only a single value within a selection,
  91  *    that value is returned over and over again.
  92  * This transparently provides the semantics of using a single selection/vector
  93  * to compute angles against multiple selections/vectors as described in the
  94  * tool help text.
  95  *
  96  * This class isn't perferctly self-contained and requires the caller to know
  97  * some of the internals to use it properly, but it serves its purpose for this
  98  * single analysis tool by simplifying the loops.
  99  * Some methods have also been tailored to allow the caller to use it a bit
 100  * more easily.
 101  *
 102  * \ingroup module_trajectoryanalysis
 103  */
 104 class AnglePositionIterator
 105 {
 106     public:
 107         /*! \brief
 108          * Creates an iterator to loop over input selection positions.
 109          *
 110          * \param[in] selections       List of selections.
 111          * \param[in] posCountPerValue Number of selection positions that
 112          *     constitute a single value for the iteration.
 113          *
 114          * If \p selections is empty, and/or \p posCountPerValue is zero, the
 115          * iterator can still be advanced and hasValue()/hasSingleValue()
 116          * called, but values cannot be accessed.
 117          */
 118         AnglePositionIterator(const SelectionList &selections,
 119                               int                  posCountPerValue)
 120             : selections_(selections), posCountPerValue_(posCountPerValue),
 121               currentSelection_(0), nextPosition_(0)
 122         {
 123         }
 124
 125         //! Advances the iterator to the next group of angles.
 126         void nextGroup()
 127         {
 128             if (selections_.size() > 1)
 129             {
 130                 ++currentSelection_;
 131             }
 132             nextPosition_ = 0;
 133         }
 134         //! Advances the iterator to the next angle in the current group.
 135         void nextValue()
 136         {
 137             if (!hasSingleValue())
 138             {
 139                 nextPosition_ += posCountPerValue_;
 140             }
 141         }
 142
 143         /*! \brief
 144          * Returns whether this iterator represents any values.
 145          *
 146          * If the return value is `false`, only nextGroup(), nextValue() and
 147          * hasSingleValue() are allowed to be called.
 148          */
 149         bool hasValue() const
 150         {
 151             return !selections_.empty();
 152         }
 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
 165         {
 166             return currentSelection().isDynamic();
 167         }
 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_ < static_cast<int>(selections_.size()),
 203                        "Accessing an invalid selection");
 204             return selections_[currentSelection_];
 205         }
 206         //! Returns the `i`th position for the current value.
 207         SelectionPosition currentPosition(int i) const
 208         {
 209             return currentSelection().position(nextPosition_ + i);
 210         }
 211         /*! \brief
 212          * Extracts all coordinates corresponding to the current value.
 213          *
 214          * \param[out] x  Array to which the positions are extracted.
 215          *
 216          * \p x should contain at minimum the number of positions per value
 217          * passed to the constructor.
 218          */
 219         void getCurrentPositions(rvec x[]) const
 220         {
 221             GMX_ASSERT(posCountPerValue_ > 0,
 222                        "Accessing positions for an invalid angle type");
 223             GMX_ASSERT(nextPosition_ + posCountPerValue_ <= currentSelection().posCount(),
 224                        "Accessing an invalid position");
 225             for (int i = 0; i < posCountPerValue_; ++i)
 226             {
 227                 copy_rvec(currentPosition(i).x(), x[i]);
 228             }
 229         }
 230
 231     private:
 232         const SelectionList &selections_;
 233         const int            posCountPerValue_;
 234         int                  currentSelection_;
 235         int                  nextPosition_;
 236
 237         GMX_DISALLOW_COPY_AND_ASSIGN(AnglePositionIterator);
 238 };
 239
 240 /********************************************************************
 241  * Actual analysis module
 242  */
 243
 244 //! How to interpret the selections in -group1.
 245 enum Group1Type
 246 {
 247     Group1Type_Angle,
 248     Group1Type_Dihedral,
 249     Group1Type_Vector,
 250     Group1Type_Plane
 251 };
 252 //! How to interpret the selections in -group2.
 253 enum Group2Type
 254 {
 255     Group2Type_None,
 256     Group2Type_Vector,
 257     Group2Type_Plane,
 258     Group2Type_TimeZero,
 259     Group2Type_Z,
 260     Group2Type_SphereNormal
 261 };
 262 //! String values corresponding to Group1Type.
 263 const char *const cGroup1TypeEnum[] =
 264 { "angle", "dihedral", "vector", "plane" };
 265 //! String values corresponding to Group2Type.
 266 const char *const cGroup2TypeEnum[] =
 267 { "none", "vector", "plane", "t0", "z", "sphnorm" };
 268
 269 class Angle : public TrajectoryAnalysisModule
 270 {
 271     public:
 272         Angle();
 273
 274         virtual void initOptions(IOptionsContainer          *options,
 275                                  TrajectoryAnalysisSettings *settings);
 276         virtual void optionsFinished(TrajectoryAnalysisSettings *settings);
 277         virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
 278                                   const TopologyInformation        &top);
 279
 280         virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 281                                   TrajectoryAnalysisModuleData *pdata);
 282
 283         virtual void finishAnalysis(int nframes);
 284         virtual void writeOutput();
 285
 286     private:
 287         void initFromSelections(const SelectionList &sel1,
 288                                 const SelectionList &sel2);
 289         void checkSelections(const SelectionList &sel1,
 290                              const SelectionList &sel2) const;
 291
 292         SelectionList                            sel1_;
 293         SelectionList                            sel2_;
 294         SelectionOptionInfo                     *sel1info_;
 295         SelectionOptionInfo                     *sel2info_;
 296         std::string                              fnAverage_;
 297         std::string                              fnAll_;
 298         std::string                              fnHistogram_;
 299
 300         Group1Type                               g1type_;
 301         Group2Type                               g2type_;
 302         double                                   binWidth_;
 303
 304         AnalysisData                             angles_;
 305         AnalysisDataFrameAverageModulePointer    averageModule_;
 306         AnalysisDataSimpleHistogramModulePointer histogramModule_;
 307
 308         std::vector<int>                         angleCount_;
 309         int                                      natoms1_;
 310         int                                      natoms2_;
 311         std::vector<std::vector<RVec> >          vt0_;
 312
 313         // Copy and assign disallowed by base.
 314 };
 315
 316 Angle::Angle()
 317     : sel1info_(nullptr), sel2info_(nullptr),
 318       g1type_(Group1Type_Angle), g2type_(Group2Type_None),
 319       binWidth_(1.0), natoms1_(0), natoms2_(0)
 320 {
 321     averageModule_.reset(new AnalysisDataFrameAverageModule());
 322     angles_.addModule(averageModule_);
 323     histogramModule_.reset(new AnalysisDataSimpleHistogramModule());
 324     angles_.addModule(histogramModule_);
 325
 326     registerAnalysisDataset(&angles_, "angle");
 327     registerBasicDataset(averageModule_.get(), "average");
 328     registerBasicDataset(&histogramModule_->averager(), "histogram");
 329 }
 330
 331
 332 void
 333 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").filetype(eftPlot).outputFile()
 387                            .store(&fnAverage_).defaultBasename("angaver")
 388                            .description("Average angles as a function of time"));
 389     options->addOption(FileNameOption("oall").filetype(eftPlot).outputFile()
 390                            .store(&fnAll_).defaultBasename("angles")
 391                            .description("All angles as a function of time"));
 392     options->addOption(FileNameOption("oh").filetype(eftPlot).outputFile()
 393                            .store(&fnHistogram_).defaultBasename("anghist")
 394                            .description("Histogram of the angles"));
 395
 396     options->addOption(EnumOption<Group1Type>("g1").enumValue(cGroup1TypeEnum)
 397                            .store(&g1type_)
 398                            .description("Type of analysis/first vector group"));
 399     options->addOption(EnumOption<Group2Type>("g2").enumValue(cGroup2TypeEnum)
 400                            .store(&g2type_)
 401                            .description("Type of second vector group"));
 402     options->addOption(DoubleOption("binw").store(&binWidth_)
 403                            .description("Binwidth for -oh in degrees"));
 404
 405     sel1info_ = options->addOption(SelectionOption("group1")
 406                                        .required().dynamicMask().storeVector(&sel1_)
 407                                        .multiValue()
 408                                        .description("First analysis/vector selection"));
 409     sel2info_ = options->addOption(SelectionOption("group2")
 410                                        .dynamicMask().storeVector(&sel2_)
 411                                        .multiValue()
 412                                        .description("Second analysis/vector selection"));
 413 }
 414
 415
 416 void
 417 Angle::optionsFinished(TrajectoryAnalysisSettings * /* settings */)
 418 {
 419     const bool bSingle = (g1type_ == Group1Type_Angle || g1type_ == Group1Type_Dihedral);
 420
 421     if (bSingle && g2type_ != Group2Type_None)
 422     {
 423         GMX_THROW(InconsistentInputError("Cannot use a second group (-g2) with "
 424                                          "-g1 angle or dihedral"));
 425     }
 426     if (bSingle && sel2info_->isSet())
 427     {
 428         GMX_THROW(InconsistentInputError("Cannot provide a second selection "
 429                                          "(-group2) with -g1 angle or dihedral"));
 430     }
 431     if (!bSingle && g2type_ == Group2Type_None)
 432     {
 433         GMX_THROW(InconsistentInputError("Should specify a second group (-g2) "
 434                                          "if the first group is not an angle or a dihedral"));
 435     }
 436
 437     // Set up the number of positions per angle.
 438     switch (g1type_)
 439     {
 440         case Group1Type_Angle:    natoms1_ = 3; break;
 441         case Group1Type_Dihedral: natoms1_ = 4; break;
 442         case Group1Type_Vector:   natoms1_ = 2; break;
 443         case Group1Type_Plane:    natoms1_ = 3; break;
 444         default:
 445             GMX_THROW(InternalError("invalid -g1 value"));
 446     }
 447     switch (g2type_)
 448     {
 449         case Group2Type_None:         natoms2_ = 0; break;
 450         case Group2Type_Vector:       natoms2_ = 2; break;
 451         case Group2Type_Plane:        natoms2_ = 3; break;
 452         case Group2Type_TimeZero:     natoms2_ = 0; break;
 453         case Group2Type_Z:            natoms2_ = 0; break;
 454         case Group2Type_SphereNormal: natoms2_ = 1; break;
 455         default:
 456             GMX_THROW(InternalError("invalid -g2 value"));
 457     }
 458     if (natoms2_ == 0 && sel2info_->isSet())
 459     {
 460         GMX_THROW(InconsistentInputError("Cannot provide a second selection (-group2) with -g2 t0 or z"));
 461     }
 462     // TODO: If bSingle is not set, the second selection option should be
 463     // required.
 464 }
 465
 466
 467 void
 468 Angle::initFromSelections(const SelectionList &sel1,
 469                           const SelectionList &sel2)
 470 {
 471     const int  angleGroups         = std::max(sel1.size(), sel2.size());
 472     const bool bHasSecondSelection = natoms2_ > 0;
 473
 474     if (bHasSecondSelection && sel1.size() != sel2.size()
 475         && std::min(sel1.size(), sel2.size()) != 1)
 476     {
 477         GMX_THROW(InconsistentInputError(
 478                           "-group1 and -group2 should specify the same number of selections"));
 479     }
 480
 481     AnglePositionIterator iter1(sel1, natoms1_);
 482     AnglePositionIterator iter2(sel2, natoms2_);
 483     for (int g = 0; g < angleGroups; ++g, iter1.nextGroup(), iter2.nextGroup())
 484     {
 485         const int posCount1 = iter1.currentSelection().posCount();
 486         if (natoms1_ > 1 && posCount1 % natoms1_ != 0)
 487         {
 488             GMX_THROW(InconsistentInputError(formatString(
 489                                                      "Number of positions in selection %d in the first group not divisible by %d",
 490                                                      static_cast<int>(g + 1), natoms1_)));
 491         }
 492         const int angleCount1 = posCount1 / natoms1_;
 493         int       angleCount  = angleCount1;
 494
 495         if (bHasSecondSelection)
 496         {
 497             const int posCount2 = iter2.currentSelection().posCount();
 498             if (natoms2_ > 1 && posCount2 % natoms2_ != 0)
 499             {
 500                 GMX_THROW(InconsistentInputError(formatString(
 501                                                          "Number of positions in selection %d in the second group not divisible by %d",
 502                                                          static_cast<int>(g + 1), natoms2_)));
 503             }
 504             if (g2type_ == Group2Type_SphereNormal && posCount2 != 1)
 505             {
 506                 GMX_THROW(InconsistentInputError(
 507                                   "The second group should contain a single position with -g2 sphnorm"));
 508             }
 509
 510             const int angleCount2 = posCount2 / natoms2_;
 511             angleCount = std::max(angleCount1, angleCount2);
 512             if (angleCount1 != angleCount2
 513                 && std::min(angleCount1, angleCount2) != 1)
 514             {
 515                 GMX_THROW(InconsistentInputError(
 516                                   "Number of vectors defined by the two groups are not the same"));
 517             }
 518         }
 519         angleCount_.push_back(angleCount);
 520     }
 521 }
 522
 523
 524 void
 525 Angle::checkSelections(const SelectionList &sel1,
 526                        const SelectionList &sel2) const
 527 {
 528     AnglePositionIterator iter1(sel1, natoms1_);
 529     AnglePositionIterator iter2(sel2, natoms2_);
 530     for (size_t g = 0; g < angleCount_.size(); ++g, iter1.nextGroup(), iter2.nextGroup())
 531     {
 532         if (iter1.isDynamic() || (iter2.hasValue() && iter2.isDynamic()))
 533         {
 534             for (int n = 0; n < angleCount_[g]; ++n, iter1.nextValue(), iter2.nextValue())
 535             {
 536                 bool bOk = true;
 537                 if (!iter1.allValuesConsistentlySelected())
 538                 {
 539                     bOk = false;
 540                 }
 541                 if (!iter2.allValuesConsistentlySelected())
 542                 {
 543                     bOk = false;
 544                 }
 545                 if (angleCount_[g] > 1)
 546                 {
 547                     if (iter1.hasSingleValue() && !iter1.currentValuesSelected())
 548                     {
 549                         bOk = false;
 550                     }
 551                     if (iter2.hasSingleValue() && !iter2.currentValuesSelected())
 552                     {
 553                         bOk = false;
 554                     }
 555                 }
 556                 if (iter2.hasValue()
 557                     && (angleCount_[g] == 1
 558                         || (!iter1.hasSingleValue() && !iter2.hasSingleValue()))
 559                     && iter1.currentValuesSelected() != iter2.currentValuesSelected())
 560                 {
 561                     bOk = false;
 562                 }
 563                 if (!bOk)
 564                 {
 565                     std::string message =
 566                         formatString("Dynamic selection %d does not select "
 567                                      "a consistent set of angles over the frames",
 568                                      static_cast<int>(g + 1));
 569                     GMX_THROW(InconsistentInputError(message));
 570                 }
 571             }
 572         }
 573     }
 574 }
 575
 576
 577 void
 578 Angle::initAnalysis(const TrajectoryAnalysisSettings &settings,
 579                     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)
 591                                .binWidth(binWidth_).includeAll());
 592
 593     if (g2type_ == Group2Type_TimeZero)
 594     {
 595         vt0_.resize(sel1_.size());
 596         for (size_t g = 0; g < sel1_.size(); ++g)
 597         {
 598             vt0_[g].resize(sel1_[g].posCount() / natoms1_);
 599         }
 600     }
 601
 602     if (!fnAverage_.empty())
 603     {
 604         AnalysisDataPlotModulePointer plotm(
 605                 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(
 622                 new AnalysisDataPlotModule(settings.plotSettings()));
 623         plotm->setFileName(fnAll_);
 624         plotm->setTitle("Angle");
 625         plotm->setXAxisIsTime();
 626         plotm->setYLabel("Angle (degrees)");
 627         // TODO: Add legends? (there can be a massive amount of columns)
 628         angles_.addModule(plotm);
 629     }
 630
 631     if (!fnHistogram_.empty())
 632     {
 633         AnalysisDataPlotModulePointer plotm(
 634                 new AnalysisDataPlotModule(settings.plotSettings()));
 635         plotm->setFileName(fnHistogram_);
 636         plotm->setTitle("Angle histogram");
 637         plotm->setXLabel("Angle (degrees)");
 638         plotm->setYLabel("Probability");
 639         // TODO: Consider adding information about the second selection,
 640         // and/or a subtitle describing what kind of angle this is.
 641         for (size_t g = 0; g < sel1_.size(); ++g)
 642         {
 643             plotm->appendLegend(sel1_[g].name());
 644         }
 645         histogramModule_->averager().addModule(plotm);
 646     }
 647 }
 648
 649
 650 //! Helper method to calculate a vector from two or three positions.
 651 static void
 652 calc_vec(int natoms, rvec x[], t_pbc *pbc, rvec xout, rvec cout)
 653 {
 654     switch (natoms)
 655     {
 656         case 2:
 657             if (pbc)
 658             {
 659                 pbc_dx(pbc, x[1], x[0], xout);
 660             }
 661             else
 662             {
 663                 rvec_sub(x[1], x[0], xout);
 664             }
 665             svmul(0.5, xout, cout);
 666             rvec_add(x[0], cout, cout);
 667             break;
 668         case 3:
 669         {
 670             rvec v1, v2;
 671             if (pbc)
 672             {
 673                 pbc_dx(pbc, x[1], x[0], v1);
 674                 pbc_dx(pbc, x[2], x[0], v2);
 675             }
 676             else
 677             {
 678                 rvec_sub(x[1], x[0], v1);
 679                 rvec_sub(x[2], x[0], v2);
 680             }
 681             cprod(v1, v2, xout);
 682             rvec_add(x[0], x[1], cout);
 683             rvec_add(cout, x[2], cout);
 684             svmul(1.0/3.0, cout, cout);
 685             break;
 686         }
 687         default:
 688             GMX_RELEASE_ASSERT(false, "Incorrectly initialized number of atoms");
 689     }
 690 }
 691
 692
 693 void
 694 Angle::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 695                     TrajectoryAnalysisModuleData *pdata)
 696 {
 697     AnalysisDataHandle       dh   = pdata->dataHandle(angles_);
 698     const SelectionList     &sel1 = pdata->parallelSelections(sel1_);
 699     const SelectionList     &sel2 = pdata->parallelSelections(sel2_);
 700
 701     checkSelections(sel1, sel2);
 702
 703     dh.startFrame(frnr, fr.time);
 704
 705     AnglePositionIterator iter1(sel1, natoms1_);
 706     AnglePositionIterator iter2(sel2, natoms2_);
 707     for (size_t g = 0; g < angleCount_.size(); ++g, iter1.nextGroup(), iter2.nextGroup())
 708     {
 709         rvec  v1, v2;
 710         rvec  c1, c2;
 711
 712         // v2 & c2 are conditionally set in the switch statement below, and conditionally
 713         // used in a different switch statement later. Apparently the clang static analyzer
 714         // thinks there are cases where they can be used uninitialzed (which I cannot find),
 715         // but to avoid trouble if we ever change just one of the switch statements it
 716         // makes sense to clear them outside the first switch.
 717
 718         clear_rvec(v2);
 719         clear_rvec(c2);
 720
 721         switch (g2type_)
 722         {
 723             case Group2Type_Z:
 724                 v2[ZZ] = 1.0;
 725                 break;
 726             case Group2Type_SphereNormal:
 727                 copy_rvec(sel2_[g].position(0).x(), c2);
 728                 break;
 729             default:
 730                 // do nothing
 731                 break;
 732         }
 733
 734         dh.selectDataSet(g);
 735         for (int n = 0; n < angleCount_[g]; ++n, iter1.nextValue(), iter2.nextValue())
 736         {
 737             rvec x[4];
 738             // x[] will be assigned below based on the number of atoms used to initialize iter1,
 739             // which in turn should correspond perfectly to g1type_ (which determines how many we read),
 740             // but unsurprisingly the static analyzer chokes a bit on that.
 741             clear_rvecs(4, x);
 742
 743             real angle;
 744             // checkSelections() ensures that this reflects all the involved
 745             // positions.
 746             const bool bPresent =
 747                 iter1.currentValuesSelected() && iter2.currentValuesSelected();
 748             iter1.getCurrentPositions(x);
 749             switch (g1type_)
 750             {
 751                 case Group1Type_Angle:
 752                     if (pbc)
 753                     {
 754                         pbc_dx(pbc, x[0], x[1], v1);
 755                         pbc_dx(pbc, x[2], x[1], v2);
 756                     }
 757                     else
 758                     {
 759                         rvec_sub(x[0], x[1], v1);
 760                         rvec_sub(x[2], x[1], v2);
 761                     }
 762                     angle = gmx_angle(v1, v2);
 763                     break;
 764                 case Group1Type_Dihedral:
 765                 {
 766                     rvec dx[3];
 767                     if (pbc)
 768                     {
 769                         // cppcheck-suppress uninitvar
 770                         pbc_dx(pbc, x[0], x[1], dx[0]);
 771                         // cppcheck-suppress uninitvar
 772                         pbc_dx(pbc, x[2], x[1], dx[1]);
 773                         // cppcheck-suppress uninitvar
 774                         pbc_dx(pbc, x[2], x[3], dx[2]);
 775                     }
 776                     else
 777                     {
 778                         // cppcheck-suppress uninitvar
 779                         rvec_sub(x[0], x[1], dx[0]);
 780                         // cppcheck-suppress uninitvar
 781                         rvec_sub(x[2], x[1], dx[1]);
 782                         // cppcheck-suppress uninitvar
 783                         rvec_sub(x[2], x[3], dx[2]);
 784                     }
 785                     cprod(dx[0], dx[1], v1);
 786                     cprod(dx[1], dx[2], v2);
 787                     angle = gmx_angle(v1, v2);
 788                     real ipr = iprod(dx[0], v2);
 789                     if (ipr < 0)
 790                     {
 791                         angle = -angle;
 792                     }
 793                     break;
 794                 }
 795                 case Group1Type_Vector:
 796                 case Group1Type_Plane:
 797                     calc_vec(natoms1_, x, pbc, v1, c1);
 798                     switch (g2type_)
 799                     {
 800                         case Group2Type_Vector:
 801                         case Group2Type_Plane:
 802                             iter2.getCurrentPositions(x);
 803                             calc_vec(natoms2_, x, pbc, v2, c2);
 804                             break;
 805                         case Group2Type_TimeZero:
 806                             // FIXME: This is not parallelizable.
 807                             if (frnr == 0)
 808                             {
 809                                 copy_rvec(v1, vt0_[g][n]);
 810                             }
 811                             copy_rvec(vt0_[g][n], v2);
 812                             break;
 813                         case Group2Type_Z:
 814                             c1[XX] = c1[YY] = 0.0;
 815                             break;
 816                         case Group2Type_SphereNormal:
 817                             if (pbc)
 818                             {
 819                                 pbc_dx(pbc, c1, c2, v2);
 820                             }
 821                             else
 822                             {
 823                                 rvec_sub(c1, c2, v2);
 824                             }
 825                             break;
 826                         default:
 827                             GMX_THROW(InternalError("invalid -g2 value"));
 828                     }
 829                     angle = gmx_angle(v1, v2);
 830                     break;
 831                 default:
 832                     GMX_THROW(InternalError("invalid -g1 value"));
 833             }
 834             dh.setPoint(n, angle * RAD2DEG, bPresent);
 835         }
 836     }
 837     dh.finishFrame();
 838 }
 839
 840
 841 void
 842 Angle::finishAnalysis(int /*nframes*/)
 843 {
 844     AbstractAverageHistogram &averageHistogram = histogramModule_->averager();
 845     averageHistogram.normalizeProbability();
 846     averageHistogram.done();
 847 }
 848
 849
 850 void
 851 Angle::writeOutput()
 852 {
 853 }
 854
 855 }       // namespace
 856
 857 const char AngleInfo::name[]             = "gangle";
 858 const char AngleInfo::shortDescription[] =
 859     "Calculate angles";
 860
 861 TrajectoryAnalysisModulePointer AngleInfo::create()
 862 {
 863     return TrajectoryAnalysisModulePointer(new Angle);
 864 }
 865
 866 } // namespace analysismodules
 867
 868 } // namespace gmx