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