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