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