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