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  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   6  * others, as listed in the AUTHORS file in the top-level source
   7  * 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 "gromacs/legacyheaders/pbc.h"
  45 #include "gromacs/legacyheaders/vec.h"
  46
  47 #include "gromacs/analysisdata/analysisdata.h"
  48 #include "gromacs/analysisdata/modules/plot.h"
  49 #include "gromacs/options/basicoptions.h"
  50 #include "gromacs/options/filenameoption.h"
  51 #include "gromacs/options/options.h"
  52 #include "gromacs/selection/selection.h"
  53 #include "gromacs/selection/selectionoption.h"
  54 #include "gromacs/trajectoryanalysis/analysissettings.h"
  55 #include "gromacs/utility/exceptions.h"
  56 #include "gromacs/utility/gmxassert.h"
  57 #include "gromacs/utility/stringutil.h"
  58
  59 namespace gmx
  60 {
  61
  62 namespace analysismodules
  63 {
  64
  65 const char Angle::name[]             = "gangle";
  66 const char Angle::shortDescription[] =
  67     "Calculate angles";
  68
  69 Angle::Angle()
  70     : TrajectoryAnalysisModule(name, shortDescription),
  71       sel1info_(NULL), sel2info_(NULL), binWidth_(1.0), natoms1_(0), natoms2_(0)
  72 {
  73     averageModule_.reset(new AnalysisDataFrameAverageModule());
  74     angles_.addModule(averageModule_);
  75     histogramModule_.reset(new AnalysisDataSimpleHistogramModule());
  76     angles_.addModule(histogramModule_);
  77
  78     registerAnalysisDataset(&angles_, "angle");
  79     registerBasicDataset(averageModule_.get(), "average");
  80     registerBasicDataset(&histogramModule_->averager(), "histogram");
  81 }
  82
  83
  84 Angle::~Angle()
  85 {
  86     for (size_t g = 0; g < vt0_.size(); ++g)
  87     {
  88         delete [] vt0_[g];
  89     }
  90 }
  91
  92
  93 void
  94 Angle::initOptions(Options *options, TrajectoryAnalysisSettings * /*settings*/)
  95 {
  96     static const char *const desc[] = {
  97         "g_angle computes different types of angles between vectors.",
  98         "It supports both vectors defined by two positions and normals of",
  99         "planes defined by three positions.",
 100         "The z axis or the local normal of a sphere can also be used as",
 101         "one of the vectors.",
 102         "There are also convenience options 'angle' and 'dihedral' for",
 103         "calculating bond angles and dihedrals defined by three/four",
 104         "positions.[PAR]",
 105         "The type of the angle is specified with [TT]-g1[tt] and [TT]-g2[tt].",
 106         "If [TT]-g1[tt] is [TT]angle[tt] or [TT]dihedral[tt], [TT]-g2[tt]",
 107         "should not be specified.",
 108         "In this case, [TT]-group1[tt] should specify one or more selections,",
 109         "and each should contain triplets or quartets of positions that define",
 110         "the angles to be calculated.[PAR]",
 111         "If [TT]-g1[tt] is [TT]vector[tt] or [TT]plane[tt], [TT]-group1[tt]",
 112         "should specify selections that contain either pairs ([TT]vector[tt])",
 113         "or triplets ([TT]plane[tt]) of positions. For vectors, the positions",
 114         "set the endpoints of the vector, and for planes, the three positions",
 115         "are used to calculate the normal of the plane. In both cases,",
 116         "[TT]-g2[tt] specifies the other vector to use (see below).[PAR]",
 117         "With [TT]-g2 vector[tt] or [TT]-g2 plane[tt], [TT]-group2[tt] should",
 118         "specify another set of vectors. [TT]-group1[tt] and [TT]-group2[tt]",
 119         "should specify the same number of selections, and for each selection",
 120         "in [TT]-group1[tt], the corresponding selection in [TT]-group2[tt]",
 121         "should specify the same number of vectors.[PAR]",
 122         "With [TT]-g2 sphnorm[tt], each selection in [TT]-group2[tt] should",
 123         "specify a single position that is the center of the sphere.",
 124         "The second vector is calculated as the vector from the center to the",
 125         "midpoint of the positions specified by [TT]-group1[tt].[PAR]",
 126         "With [TT]-g2 z[tt], [TT]-group2[tt] is not necessary, and angles",
 127         "between the first vectors and the positive Z axis are calculated.[PAR]",
 128         "With [TT]-g2 t0[tt], [TT]-group2[tt] is not necessary, and angles",
 129         "are calculated from the vectors as they are in the first frame.[PAR]",
 130         "There are three options for output:",
 131         "[TT]-oav[tt] writes an xvgr file with the time and the average angle",
 132         "for each frame.",
 133         "[TT]-oall[tt] writes all the individual angles.",
 134         "[TT]-oh[tt] writes a histogram of the angles. The bin width can be",
 135         "set with [TT]-binw[tt].",
 136         "For [TT]-oav[tt] and [TT]-oh[tt], separate average/histogram is",
 137         "computed for each selection in [TT]-group1[tt]."
 138     };
 139     static const char *const cGroup1TypeEnum[] =
 140     { "angle", "dihedral", "vector", "plane" };
 141     static const char *const cGroup2TypeEnum[] =
 142     { "none", "vector", "plane", "t0", "z", "sphnorm" };
 143
 144     options->setDescription(concatenateStrings(desc));
 145
 146     options->addOption(FileNameOption("oav").filetype(eftPlot).outputFile()
 147                            .store(&fnAverage_).defaultBasename("angaver")
 148                            .description("Average angles as a function of time"));
 149     options->addOption(FileNameOption("oall").filetype(eftPlot).outputFile()
 150                            .store(&fnAll_).defaultBasename("angles")
 151                            .description("All angles as a function of time"));
 152     options->addOption(FileNameOption("oh").filetype(eftPlot).outputFile()
 153                            .store(&fnHistogram_).defaultBasename("anghist")
 154                            .description("Histogram of the angles"));
 155
 156     options->addOption(StringOption("g1").enumValue(cGroup1TypeEnum)
 157                            .defaultEnumIndex(0).store(&g1type_)
 158                            .description("Type of analysis/first vector group"));
 159     options->addOption(StringOption("g2").enumValue(cGroup2TypeEnum)
 160                            .defaultEnumIndex(0).store(&g2type_)
 161                            .description("Type of second vector group"));
 162     options->addOption(DoubleOption("binw").store(&binWidth_)
 163                            .description("Binwidth for -oh in degrees"));
 164
 165     sel1info_ = options->addOption(SelectionOption("group1")
 166                                        .required().dynamicMask().storeVector(&sel1_)
 167                                        .multiValue()
 168                                        .description("First analysis/vector selection"));
 169     sel2info_ = options->addOption(SelectionOption("group2")
 170                                        .dynamicMask().storeVector(&sel2_)
 171                                        .multiValue()
 172                                        .description("Second analysis/vector selection"));
 173 }
 174
 175
 176 void
 177 Angle::optionsFinished(Options *options, TrajectoryAnalysisSettings *settings)
 178 {
 179     bool bSingle = (g1type_[0] == 'a' || g1type_[0] == 'd');
 180
 181     if (bSingle && g2type_[0] != 'n')
 182     {
 183         GMX_THROW(InconsistentInputError("Cannot use a second group (-g2) with "
 184                                          "-g1 angle or dihedral"));
 185     }
 186     if (bSingle && options->isSet("group2"))
 187     {
 188         GMX_THROW(InconsistentInputError("Cannot provide a second selection "
 189                                          "(-group2) with -g1 angle or dihedral"));
 190     }
 191     if (!bSingle && g2type_[0] == 'n')
 192     {
 193         GMX_THROW(InconsistentInputError("Should specify a second group (-g2) "
 194                                          "if the first group is not an angle or a dihedral"));
 195     }
 196
 197     // Set up the number of positions per angle.
 198     switch (g1type_[0])
 199     {
 200         case 'a': natoms1_ = 3; break;
 201         case 'd': natoms1_ = 4; break;
 202         case 'v': natoms1_ = 2; break;
 203         case 'p': natoms1_ = 3; break;
 204         default:
 205             GMX_THROW(InternalError("invalid -g1 value"));
 206     }
 207     switch (g2type_[0])
 208     {
 209         case 'n': natoms2_ = 0; break;
 210         case 'v': natoms2_ = 2; break;
 211         case 'p': natoms2_ = 3; break;
 212         case 't': natoms2_ = 0; break;
 213         case 'z': natoms2_ = 0; break;
 214         case 's': natoms2_ = 1; break;
 215         default:
 216             GMX_THROW(InternalError("invalid -g2 value"));
 217     }
 218     if (natoms2_ == 0 && options->isSet("group2"))
 219     {
 220         GMX_THROW(InconsistentInputError("Cannot provide a second selection (-group2) with -g2 t0 or z"));
 221     }
 222 }
 223
 224
 225 void
 226 Angle::checkSelections(const SelectionList &sel1,
 227                        const SelectionList &sel2) const
 228 {
 229     if (natoms2_ > 0 && sel1.size() != sel2.size())
 230     {
 231         GMX_THROW(InconsistentInputError(
 232                           "-group1 and -group2 should specify the same number of selections"));
 233     }
 234
 235     for (size_t g = 0; g < sel1.size(); ++g)
 236     {
 237         const int na1 = sel1[g].posCount();
 238         const int na2 = (natoms2_ > 0) ? sel2[g].posCount() : 0;
 239         if (natoms1_ > 1 && na1 % natoms1_ != 0)
 240         {
 241             GMX_THROW(InconsistentInputError(formatString(
 242                                                      "Number of positions in selection %d in the first group not divisible by %d",
 243                                                      static_cast<int>(g + 1), natoms1_)));
 244         }
 245         if (natoms2_ > 1 && na2 % natoms2_ != 0)
 246         {
 247             GMX_THROW(InconsistentInputError(formatString(
 248                                                      "Number of positions in selection %d in the second group not divisible by %d",
 249                                                      static_cast<int>(g + 1), natoms2_)));
 250         }
 251         if (natoms1_ > 0 && natoms2_ > 1 && na1 / natoms1_ != na2 / natoms2_)
 252         {
 253             GMX_THROW(InconsistentInputError(
 254                               "Number of vectors defined by the two groups are not the same"));
 255         }
 256         if (g2type_[0] == 's' && sel2[g].posCount() != 1)
 257         {
 258             GMX_THROW(InconsistentInputError(
 259                               "The second group should contain a single position with -g2 sphnorm"));
 260         }
 261         if (sel1[g].isDynamic() || (natoms2_ > 0 && sel2[g].isDynamic()))
 262         {
 263             for (int i = 0, j = 0; i < na1; i += natoms1_, j += natoms2_)
 264             {
 265                 const bool bSelected = sel1[g].position(i).selected();
 266                 bool       bOk       = true;
 267                 for (int k = 1; k < natoms1_ && bOk; ++k)
 268                 {
 269                     bOk = (sel1[g].position(i+k).selected() == bSelected);
 270                 }
 271                 for (int k = 1; k < natoms2_ && bOk; ++k)
 272                 {
 273                     bOk = (sel2[g].position(j+k).selected() == bSelected);
 274                 }
 275                 if (!bOk)
 276                 {
 277                     std::string message =
 278                         formatString("Dynamic selection %d does not select "
 279                                      "a consistent set of angles over the frames",
 280                                      static_cast<int>(g + 1));
 281                     GMX_THROW(InconsistentInputError(message));
 282                 }
 283             }
 284         }
 285     }
 286 }
 287
 288
 289 void
 290 Angle::initAnalysis(const TrajectoryAnalysisSettings &settings,
 291                     const TopologyInformation        &top)
 292 {
 293     checkSelections(sel1_, sel2_);
 294
 295     // checkSelections() ensures that both selection lists have the same size.
 296     angles_.setDataSetCount(sel1_.size());
 297     for (size_t i = 0; i < sel1_.size(); ++i)
 298     {
 299         angles_.setColumnCount(i, sel1_[i].posCount() / natoms1_);
 300     }
 301     double histogramMin = (g1type_ == "dihedral" ? -180.0 : 0);
 302     histogramModule_->init(histogramFromRange(histogramMin, 180.0)
 303                                .binWidth(binWidth_).includeAll());
 304
 305     if (g2type_ == "t0")
 306     {
 307         vt0_.resize(sel1_.size());
 308         for (size_t g = 0; g < sel1_.size(); ++g)
 309         {
 310             vt0_[g] = new rvec[sel1_[g].posCount() / natoms1_];
 311         }
 312     }
 313
 314     if (!fnAverage_.empty())
 315     {
 316         AnalysisDataPlotModulePointer plotm(
 317                 new AnalysisDataPlotModule(settings.plotSettings()));
 318         plotm->setFileName(fnAverage_);
 319         plotm->setTitle("Average angle");
 320         plotm->setXAxisIsTime();
 321         plotm->setYLabel("Angle (degrees)");
 322         // TODO: Consider adding information about the second selection,
 323         // and/or a subtitle describing what kind of angle this is.
 324         for (size_t g = 0; g < sel1_.size(); ++g)
 325         {
 326             plotm->appendLegend(sel1_[g].name());
 327         }
 328         averageModule_->addModule(plotm);
 329     }
 330
 331     if (!fnAll_.empty())
 332     {
 333         AnalysisDataPlotModulePointer plotm(
 334                 new AnalysisDataPlotModule(settings.plotSettings()));
 335         plotm->setFileName(fnAll_);
 336         plotm->setTitle("Angle");
 337         plotm->setXAxisIsTime();
 338         plotm->setYLabel("Angle (degrees)");
 339         // TODO: Add legends? (there can be a massive amount of columns)
 340         angles_.addModule(plotm);
 341     }
 342
 343     if (!fnHistogram_.empty())
 344     {
 345         AnalysisDataPlotModulePointer plotm(
 346                 new AnalysisDataPlotModule(settings.plotSettings()));
 347         plotm->setFileName(fnHistogram_);
 348         plotm->setTitle("Angle histogram");
 349         plotm->setXLabel("Angle (degrees)");
 350         plotm->setYLabel("Probability");
 351         // TODO: Consider adding information about the second selection,
 352         // and/or a subtitle describing what kind of angle this is.
 353         for (size_t g = 0; g < sel1_.size(); ++g)
 354         {
 355             plotm->appendLegend(sel1_[g].name());
 356         }
 357         histogramModule_->averager().addModule(plotm);
 358     }
 359 }
 360
 361
 362 //! Helper method to process selections into an array of coordinates.
 363 static void
 364 copy_pos(const SelectionList &sel, int natoms, int g, int first, rvec x[])
 365 {
 366     for (int k = 0; k < natoms; ++k)
 367     {
 368         copy_rvec(sel[g].position(first + k).x(), x[k]);
 369     }
 370 }
 371
 372
 373 //! Helper method to calculate a vector from two or three positions.
 374 static void
 375 calc_vec(int natoms, rvec x[], t_pbc *pbc, rvec xout, rvec cout)
 376 {
 377     switch (natoms)
 378     {
 379         case 2:
 380             if (pbc)
 381             {
 382                 pbc_dx(pbc, x[1], x[0], xout);
 383             }
 384             else
 385             {
 386                 rvec_sub(x[1], x[0], xout);
 387             }
 388             svmul(0.5, xout, cout);
 389             rvec_add(x[0], cout, cout);
 390             break;
 391         case 3:
 392         {
 393             rvec v1, v2;
 394             if (pbc)
 395             {
 396                 pbc_dx(pbc, x[1], x[0], v1);
 397                 pbc_dx(pbc, x[2], x[0], v2);
 398             }
 399             else
 400             {
 401                 rvec_sub(x[1], x[0], v1);
 402                 rvec_sub(x[2], x[0], v2);
 403             }
 404             cprod(v1, v2, xout);
 405             rvec_add(x[0], x[1], cout);
 406             rvec_add(cout, x[2], cout);
 407             svmul(1.0/3.0, cout, cout);
 408             break;
 409         }
 410         default:
 411             GMX_RELEASE_ASSERT(false, "Incorrectly initialized number of atoms");
 412     }
 413 }
 414
 415
 416 void
 417 Angle::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 418                     TrajectoryAnalysisModuleData *pdata)
 419 {
 420     AnalysisDataHandle       dh   = pdata->dataHandle(angles_);
 421     const SelectionList     &sel1 = pdata->parallelSelections(sel1_);
 422     const SelectionList     &sel2 = pdata->parallelSelections(sel2_);
 423
 424     checkSelections(sel1, sel2);
 425
 426     dh.startFrame(frnr, fr.time);
 427
 428     for (size_t g = 0; g < sel1_.size(); ++g)
 429     {
 430         rvec  v1, v2;
 431         rvec  c1, c2;
 432         switch (g2type_[0])
 433         {
 434             case 'z':
 435                 clear_rvec(v2);
 436                 v2[ZZ] = 1.0;
 437                 clear_rvec(c2);
 438                 break;
 439             case 's':
 440                 copy_rvec(sel2_[g].position(0).x(), c2);
 441                 break;
 442         }
 443         dh.selectDataSet(g);
 444         for (int i = 0, j = 0, n = 0;
 445              i < sel1[g].posCount();
 446              i += natoms1_, j += natoms2_, ++n)
 447         {
 448             rvec x[4];
 449             real angle;
 450             // checkSelections() ensures that this reflects all the involved
 451             // positions.
 452             bool bPresent = sel1[g].position(i).selected();
 453             copy_pos(sel1, natoms1_, g, i, x);
 454             switch (g1type_[0])
 455             {
 456                 case 'a':
 457                     if (pbc)
 458                     {
 459                         pbc_dx(pbc, x[0], x[1], v1);
 460                         pbc_dx(pbc, x[2], x[1], v2);
 461                     }
 462                     else
 463                     {
 464                         rvec_sub(x[0], x[1], v1);
 465                         rvec_sub(x[2], x[1], v2);
 466                     }
 467                     angle = gmx_angle(v1, v2);
 468                     break;
 469                 case 'd':
 470                 {
 471                     rvec dx[3];
 472                     if (pbc)
 473                     {
 474                         pbc_dx(pbc, x[0], x[1], dx[0]);
 475                         pbc_dx(pbc, x[2], x[1], dx[1]);
 476                         pbc_dx(pbc, x[2], x[3], dx[2]);
 477                     }
 478                     else
 479                     {
 480                         rvec_sub(x[0], x[1], dx[0]);
 481                         rvec_sub(x[2], x[1], dx[1]);
 482                         rvec_sub(x[2], x[3], dx[2]);
 483                     }
 484                     cprod(dx[0], dx[1], v1);
 485                     cprod(dx[1], dx[2], v2);
 486                     angle = gmx_angle(v1, v2);
 487                     real ipr = iprod(dx[0], v2);
 488                     if (ipr < 0)
 489                     {
 490                         angle = -angle;
 491                     }
 492                     break;
 493                 }
 494                 case 'v':
 495                 case 'p':
 496                     calc_vec(natoms1_, x, pbc, v1, c1);
 497                     switch (g2type_[0])
 498                     {
 499                         case 'v':
 500                         case 'p':
 501                             copy_pos(sel2, natoms2_, 0, j, x);
 502                             calc_vec(natoms2_, x, pbc, v2, c2);
 503                             break;
 504                         case 't':
 505                             // FIXME: This is not parallelizable.
 506                             if (frnr == 0)
 507                             {
 508                                 copy_rvec(v1, vt0_[g][n]);
 509                             }
 510                             copy_rvec(vt0_[g][n], v2);
 511                             break;
 512                         case 'z':
 513                             c1[XX] = c1[YY] = 0.0;
 514                             break;
 515                         case 's':
 516                             if (pbc)
 517                             {
 518                                 pbc_dx(pbc, c1, c2, v2);
 519                             }
 520                             else
 521                             {
 522                                 rvec_sub(c1, c2, v2);
 523                             }
 524                             break;
 525                         default:
 526                             GMX_THROW(InternalError("invalid -g2 value"));
 527                     }
 528                     angle = gmx_angle(v1, v2);
 529                     break;
 530                 default:
 531                     GMX_THROW(InternalError("invalid -g1 value"));
 532             }
 533             dh.setPoint(n, angle * RAD2DEG, bPresent);
 534         }
 535     }
 536     dh.finishFrame();
 537 }
 538
 539
 540 void
 541 Angle::finishAnalysis(int /*nframes*/)
 542 {
 543     AbstractAverageHistogram &averageHistogram = histogramModule_->averager();
 544     averageHistogram.normalizeProbability();
 545     averageHistogram.done();
 546 }
 547
 548
 549 void
 550 Angle::writeOutput()
 551 {
 552 }
 553
 554 } // namespace modules
 555
 556 } // namespace gmxana