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