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