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