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: Add legends
 323         averageModule_->addModule(plotm);
 324     }
 325
 326     if (!fnAll_.empty())
 327     {
 328         AnalysisDataPlotModulePointer plotm(
 329                 new AnalysisDataPlotModule(settings.plotSettings()));
 330         plotm->setFileName(fnAll_);
 331         plotm->setTitle("Angle");
 332         plotm->setXAxisIsTime();
 333         plotm->setYLabel("Angle (degrees)");
 334         // TODO: Add legends? (there can be a massive amount of columns)
 335         angles_.addModule(plotm);
 336     }
 337
 338     if (!fnHistogram_.empty())
 339     {
 340         AnalysisDataPlotModulePointer plotm(
 341                 new AnalysisDataPlotModule(settings.plotSettings()));
 342         plotm->setFileName(fnHistogram_);
 343         plotm->setTitle("Angle histogram");
 344         plotm->setXLabel("Angle (degrees)");
 345         plotm->setYLabel("Probability");
 346         // TODO: Add legends
 347         histogramModule_->averager().addModule(plotm);
 348     }
 349 }
 350
 351
 352 //! Helper method to process selections into an array of coordinates.
 353 static void
 354 copy_pos(const SelectionList &sel, int natoms, int g, int first, rvec x[])
 355 {
 356     for (int k = 0; k < natoms; ++k)
 357     {
 358         copy_rvec(sel[g].position(first + k).x(), x[k]);
 359     }
 360 }
 361
 362
 363 //! Helper method to calculate a vector from two or three positions.
 364 static void
 365 calc_vec(int natoms, rvec x[], t_pbc *pbc, rvec xout, rvec cout)
 366 {
 367     switch (natoms)
 368     {
 369         case 2:
 370             if (pbc)
 371             {
 372                 pbc_dx(pbc, x[1], x[0], xout);
 373             }
 374             else
 375             {
 376                 rvec_sub(x[1], x[0], xout);
 377             }
 378             svmul(0.5, xout, cout);
 379             rvec_add(x[0], cout, cout);
 380             break;
 381         case 3:
 382         {
 383             rvec v1, v2;
 384             if (pbc)
 385             {
 386                 pbc_dx(pbc, x[1], x[0], v1);
 387                 pbc_dx(pbc, x[2], x[0], v2);
 388             }
 389             else
 390             {
 391                 rvec_sub(x[1], x[0], v1);
 392                 rvec_sub(x[2], x[0], v2);
 393             }
 394             cprod(v1, v2, xout);
 395             rvec_add(x[0], x[1], cout);
 396             rvec_add(cout, x[2], cout);
 397             svmul(1.0/3.0, cout, cout);
 398             break;
 399         }
 400         default:
 401             GMX_RELEASE_ASSERT(false, "Incorrectly initialized number of atoms");
 402     }
 403 }
 404
 405
 406 void
 407 Angle::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 408                     TrajectoryAnalysisModuleData *pdata)
 409 {
 410     AnalysisDataHandle       dh   = pdata->dataHandle(angles_);
 411     const SelectionList     &sel1 = pdata->parallelSelections(sel1_);
 412     const SelectionList     &sel2 = pdata->parallelSelections(sel2_);
 413
 414     checkSelections(sel1, sel2);
 415
 416     dh.startFrame(frnr, fr.time);
 417
 418     for (size_t g = 0; g < sel1_.size(); ++g)
 419     {
 420         rvec  v1, v2;
 421         rvec  c1, c2;
 422         switch (g2type_[0])
 423         {
 424             case 'z':
 425                 clear_rvec(v2);
 426                 v2[ZZ] = 1.0;
 427                 clear_rvec(c2);
 428                 break;
 429             case 's':
 430                 copy_rvec(sel2_[g].position(0).x(), c2);
 431                 break;
 432         }
 433         dh.selectDataSet(g);
 434         for (int i = 0, j = 0, n = 0;
 435              i < sel1[g].posCount();
 436              i += natoms1_, j += natoms2_, ++n)
 437         {
 438             rvec x[4];
 439             real angle;
 440             // checkSelections() ensures that this reflects all the involved
 441             // positions.
 442             bool bPresent = sel1[g].position(i).selected();
 443             copy_pos(sel1, natoms1_, g, i, x);
 444             switch (g1type_[0])
 445             {
 446                 case 'a':
 447                     if (pbc)
 448                     {
 449                         pbc_dx(pbc, x[0], x[1], v1);
 450                         pbc_dx(pbc, x[2], x[1], v2);
 451                     }
 452                     else
 453                     {
 454                         rvec_sub(x[0], x[1], v1);
 455                         rvec_sub(x[2], x[1], v2);
 456                     }
 457                     angle = gmx_angle(v1, v2);
 458                     break;
 459                 case 'd':
 460                 {
 461                     rvec dx[3];
 462                     if (pbc)
 463                     {
 464                         pbc_dx(pbc, x[0], x[1], dx[0]);
 465                         pbc_dx(pbc, x[2], x[1], dx[1]);
 466                         pbc_dx(pbc, x[2], x[3], dx[2]);
 467                     }
 468                     else
 469                     {
 470                         rvec_sub(x[0], x[1], dx[0]);
 471                         rvec_sub(x[2], x[1], dx[1]);
 472                         rvec_sub(x[2], x[3], dx[2]);
 473                     }
 474                     cprod(dx[0], dx[1], v1);
 475                     cprod(dx[1], dx[2], v2);
 476                     angle = gmx_angle(v1, v2);
 477                     real ipr = iprod(dx[0], v2);
 478                     if (ipr < 0)
 479                     {
 480                         angle = -angle;
 481                     }
 482                     break;
 483                 }
 484                 case 'v':
 485                 case 'p':
 486                     calc_vec(natoms1_, x, pbc, v1, c1);
 487                     switch (g2type_[0])
 488                     {
 489                         case 'v':
 490                         case 'p':
 491                             copy_pos(sel2, natoms2_, 0, j, x);
 492                             calc_vec(natoms2_, x, pbc, v2, c2);
 493                             break;
 494                         case 't':
 495                             // FIXME: This is not parallelizable.
 496                             if (frnr == 0)
 497                             {
 498                                 copy_rvec(v1, vt0_[g][n]);
 499                             }
 500                             copy_rvec(vt0_[g][n], v2);
 501                             break;
 502                         case 'z':
 503                             c1[XX] = c1[YY] = 0.0;
 504                             break;
 505                         case 's':
 506                             if (pbc)
 507                             {
 508                                 pbc_dx(pbc, c1, c2, v2);
 509                             }
 510                             else
 511                             {
 512                                 rvec_sub(c1, c2, v2);
 513                             }
 514                             break;
 515                         default:
 516                             GMX_THROW(InternalError("invalid -g2 value"));
 517                     }
 518                     angle = gmx_angle(v1, v2);
 519                     break;
 520                 default:
 521                     GMX_THROW(InternalError("invalid -g1 value"));
 522             }
 523             dh.setPoint(n, angle * RAD2DEG, bPresent);
 524         }
 525     }
 526     dh.finishFrame();
 527 }
 528
 529
 530 void
 531 Angle::finishAnalysis(int /*nframes*/)
 532 {
 533     AbstractAverageHistogram &averageHistogram = histogramModule_->averager();
 534     averageHistogram.normalizeProbability();
 535     averageHistogram.done();
 536 }
 537
 538
 539 void
 540 Angle::writeOutput()
 541 {
 542 }
 543
 544 } // namespace modules
 545
 546 } // namespace gmxana