src/gromacs/trajectoryanalysis/modules/angle.cpp

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  10  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  11  * Copyright (c) 2001-2009, The GROMACS development team,
  12  * check out http://www.gromacs.org for more information.
  13
  14  * This program is free software; you can redistribute it and/or
  15  * modify it under the terms of the GNU General Public License
  16  * as published by the Free Software Foundation; either version 2
  17  * of the License, or (at your option) any later version.
  18  *
  19  * If you want to redistribute modifications, please consider that
  20  * scientific software is very special. Version control is crucial -
  21  * bugs must be traceable. We will be happy to consider code for
  22  * inclusion in the official distribution, but derived work must not
  23  * be called official GROMACS. Details are found in the README & COPYING
  24  * files - if they are missing, get the official version at www.gromacs.org.
  25  *
  26  * To help us fund GROMACS development, we humbly ask that you cite
  27  * the papers on the package - you can find them in the top README file.
  28  *
  29  * For more info, check our website at http://www.gromacs.org
  30  */
  31 /*! \internal \file
  32  * \brief
  33  * Implements gmx::analysismodules::Angle.
  34  *
  35  * \author Teemu Murtola <teemu.murtola@cbr.su.se>
  36  * \ingroup module_trajectoryanalysis
  37  */
  38 #include "angle.h"
  39
  40 #include "gromacs/legacyheaders/pbc.h"
  41 #include "gromacs/legacyheaders/vec.h"
  42
  43 #include "gromacs/analysisdata/analysisdata.h"
  44 #include "gromacs/analysisdata/modules/plot.h"
  45 #include "gromacs/options/basicoptions.h"
  46 #include "gromacs/options/filenameoption.h"
  47 #include "gromacs/options/options.h"
  48 #include "gromacs/selection/selection.h"
  49 #include "gromacs/selection/selectionoption.h"
  50 #include "gromacs/trajectoryanalysis/analysissettings.h"
  51 #include "gromacs/utility/exceptions.h"
  52 #include "gromacs/utility/gmxassert.h"
  53 #include "gromacs/utility/stringutil.h"
  54
  55 namespace gmx
  56 {
  57
  58 namespace analysismodules
  59 {
  60
  61 const char Angle::name[]             = "gangle";
  62 const char Angle::shortDescription[] =
  63     "Calculate angles";
  64
  65 Angle::Angle()
  66     : TrajectoryAnalysisModule(name, shortDescription),
  67       sel1info_(NULL), sel2info_(NULL), natoms1_(0), natoms2_(0)
  68 {
  69     averageModule_.reset(new AnalysisDataFrameAverageModule());
  70     angles_.addModule(averageModule_);
  71
  72     registerAnalysisDataset(&angles_, "angle");
  73     registerBasicDataset(averageModule_.get(), "average");
  74 }
  75
  76
  77 Angle::~Angle()
  78 {
  79     for (size_t g = 0; g < vt0_.size(); ++g)
  80     {
  81         delete [] vt0_[g];
  82     }
  83 }
  84
  85
  86 void
  87 Angle::initOptions(Options *options, TrajectoryAnalysisSettings * /*settings*/)
  88 {
  89     static const char *const desc[] = {
  90         "g_angle computes different types of angles between vectors.",
  91         "It supports both vectors defined by two positions and normals of",
  92         "planes defined by three positions.",
  93         "The z axis or the local normal of a sphere can also be used as",
  94         "one of the vectors.",
  95         "There are also convenience options 'angle' and 'dihedral' for",
  96         "calculating bond angles and dihedrals defined by three/four",
  97         "positions.[PAR]",
  98         "The type of the angle is specified with [TT]-g1[tt] and [TT]-g2[tt].",
  99         "If [TT]-g1[tt] is [TT]angle[tt] or [TT]dihedral[tt], [TT]-g2[tt]",
 100         "should not be specified.",
 101         "In this case, [TT]-group1[tt] should specify one selection,",
 102         "and it should contain triplets or quartets of positions that define",
 103         "the angles to be calculated.[PAR]",
 104         "If [TT]-g1[tt] is [TT]vector[tt] or [TT]plane[tt], [TT]-group1[tt]",
 105         "should specify a selection that has either pairs ([TT]vector[tt])",
 106         "or triplets ([TT]plane[tt]) of positions. For vectors, the positions",
 107         "set the endpoints of the vector, and for planes, the three positions",
 108         "are used to calculate the normal of the plane. In both cases,",
 109         "[TT]-g2[tt] specifies the other vector to use (see below).[PAR]",
 110         "With [TT]-g2 vector[tt] or [TT]-g2 plane[tt], [TT]-group2[tt] should",
 111         "specify another set of vectors. Both selections should specify the",
 112         "same number of vectors.[PAR]",
 113         "With [TT]-g2 sphnorm[tt], [TT]-group2[tt] should specify a single",
 114         "position that is the center of the sphere. The second vector is then",
 115         "calculated as the vector from the center to the midpoint of the",
 116         "positions specified by [TT]-group1[tt].[PAR]",
 117         "With [TT]-g2 z[tt], [TT]-group2[tt] is not necessary, and angles",
 118         "between the first vectors and the positive Z axis are calculated.[PAR]",
 119         "With [TT]-g2 t0[tt], [TT]-group2[tt] is not necessary, and angles",
 120         "are calculated from the vectors as they are in the first frame.[PAR]",
 121         "There are two options for output:",
 122         "[TT]-oav[tt] writes an xvgr file with the time and the average angle",
 123         "for each frame.",
 124         "[TT]-oall[tt] writes all the individual angles."
 125         /* TODO: Consider if the dump option is necessary and how to best
 126          * implement it.
 127            "[TT]-od[tt] can be used to dump all the individual angles,",
 128            "each on a separate line. This format is better suited for",
 129            "further processing, e.g., if angles from multiple runs are needed."
 130          */
 131     };
 132     static const char *const cGroup1TypeEnum[] =
 133     { "angle", "dihedral", "vector", "plane", NULL };
 134     static const char *const cGroup2TypeEnum[] =
 135     { "none", "vector", "plane", "t0", "z", "sphnorm", NULL };
 136
 137     options->setDescription(concatenateStrings(desc));
 138
 139     options->addOption(FileNameOption("oav").filetype(eftPlot).outputFile()
 140                            .store(&fnAverage_).defaultBasename("angaver")
 141                            .description("Average angles as a function of time"));
 142     options->addOption(FileNameOption("oall").filetype(eftPlot).outputFile()
 143                            .store(&fnAll_).defaultBasename("angles")
 144                            .description("All angles as a function of time"));
 145     // TODO: Add histogram output.
 146
 147     options->addOption(StringOption("g1").enumValue(cGroup1TypeEnum)
 148                            .defaultEnumIndex(0).store(&g1type_)
 149                            .description("Type of analysis/first vector group"));
 150     options->addOption(StringOption("g2").enumValue(cGroup2TypeEnum)
 151                            .defaultEnumIndex(0).store(&g2type_)
 152                            .description("Type of second vector group"));
 153
 154     // TODO: Allow multiple angles to be computed in one invocation.
 155     // Most of the code already supports it, but requires a solution for
 156     // Redmine issue #1010.
 157     // TODO: Consider what is the best way to support dynamic selections.
 158     // Again, most of the code already supports it, but it needs to be
 159     // considered how should -oall work, and additional checks should be added.
 160     sel1info_ = options->addOption(SelectionOption("group1")
 161                                        .required().onlyStatic().storeVector(&sel1_)
 162                                        .description("First analysis/vector selection"));
 163     sel2info_ = options->addOption(SelectionOption("group2")
 164                                        .onlyStatic().storeVector(&sel2_)
 165                                        .description("Second analysis/vector selection"));
 166 }
 167
 168
 169 void
 170 Angle::optionsFinished(Options *options, TrajectoryAnalysisSettings *settings)
 171 {
 172     bool bSingle = (g1type_[0] == 'a' || g1type_[0] == 'd');
 173
 174     if (bSingle && g2type_[0] != 'n')
 175     {
 176         GMX_THROW(InconsistentInputError("Cannot use a second group (-g2) with "
 177                                          "-g1 angle or dihedral"));
 178     }
 179     if (bSingle && options->isSet("group2"))
 180     {
 181         GMX_THROW(InconsistentInputError("Cannot provide a second selection "
 182                                          "(-group2) with -g1 angle or dihedral"));
 183     }
 184     if (!bSingle && g2type_[0] == 'n')
 185     {
 186         GMX_THROW(InconsistentInputError("Should specify a second group (-g2) "
 187                                          "if the first group is not an angle or a dihedral"));
 188     }
 189
 190     // Set up the number of positions per angle.
 191     switch (g1type_[0])
 192     {
 193         case 'a': natoms1_ = 3; break;
 194         case 'd': natoms1_ = 4; break;
 195         case 'v': natoms1_ = 2; break;
 196         case 'p': natoms1_ = 3; break;
 197         default:
 198             GMX_THROW(InternalError("invalid -g1 value"));
 199     }
 200     switch (g2type_[0])
 201     {
 202         case 'n': natoms2_ = 0; break;
 203         case 'v': natoms2_ = 2; break;
 204         case 'p': natoms2_ = 3; break;
 205         case 't': natoms2_ = 0; break;
 206         case 'z': natoms2_ = 0; break;
 207         case 's': natoms2_ = 1; break;
 208         default:
 209             GMX_THROW(InternalError("invalid -g2 value"));
 210     }
 211     if (natoms2_ == 0 && options->isSet("group2"))
 212     {
 213         GMX_THROW(InconsistentInputError("Cannot provide a second selection (-group2) with -g2 t0 or z"));
 214     }
 215 }
 216
 217
 218 void
 219 Angle::checkSelections(const SelectionList &sel1,
 220                        const SelectionList &sel2) const
 221 {
 222     if (natoms2_ > 0 && sel1.size() != sel2.size())
 223     {
 224         GMX_THROW(InconsistentInputError(
 225                           "-group1 and -group2 should specify the same number of selections"));
 226     }
 227
 228     for (size_t g = 0; g < sel1.size(); ++g)
 229     {
 230         int na1 = sel1[g].posCount();
 231         int na2 = (natoms2_ > 0) ? sel2[g].posCount() : 0;
 232         if (natoms1_ > 1 && na1 % natoms1_ != 0)
 233         {
 234             GMX_THROW(InconsistentInputError(formatString(
 235                                                      "Number of positions in selection %d in the first group not divisible by %d",
 236                                                      static_cast<int>(g + 1), natoms1_)));
 237         }
 238         if (natoms2_ > 1 && na2 % natoms2_ != 0)
 239         {
 240             GMX_THROW(InconsistentInputError(formatString(
 241                                                      "Number of positions in selection %d in the second group not divisible by %d",
 242                                                      static_cast<int>(g + 1), natoms2_)));
 243         }
 244         if (natoms1_ > 0 && natoms2_ > 1 && na1 / natoms1_ != na2 / natoms2_)
 245         {
 246             GMX_THROW(InconsistentInputError(
 247                               "Number of vectors defined by the two groups are not the same"));
 248         }
 249         if (g2type_[0] == 's' && sel2[g].posCount() != 1)
 250         {
 251             GMX_THROW(InconsistentInputError(
 252                               "The second group should contain a single position with -g2 sphnorm"));
 253         }
 254     }
 255 }
 256
 257
 258 void
 259 Angle::initAnalysis(const TrajectoryAnalysisSettings &settings,
 260                     const TopologyInformation        &top)
 261 {
 262     checkSelections(sel1_, sel2_);
 263
 264     angles_.setColumnCount(sel1_[0].posCount() / natoms1_);
 265
 266     if (g2type_ == "t0")
 267     {
 268         vt0_.resize(sel1_.size());
 269         for (size_t g = 0; g < sel1_.size(); ++g)
 270         {
 271             vt0_[g] = new rvec[sel1_[g].posCount() / natoms1_];
 272         }
 273     }
 274
 275     if (!fnAverage_.empty())
 276     {
 277         AnalysisDataPlotModulePointer plotm(
 278                 new AnalysisDataPlotModule(settings.plotSettings()));
 279         plotm->setFileName(fnAverage_);
 280         plotm->setTitle("Average angle");
 281         plotm->setXAxisIsTime();
 282         plotm->setYLabel("Angle (degrees)");
 283         // TODO: Add legends
 284         averageModule_->addModule(plotm);
 285     }
 286
 287     if (!fnAll_.empty())
 288     {
 289         AnalysisDataPlotModulePointer plotm(
 290                 new AnalysisDataPlotModule(settings.plotSettings()));
 291         plotm->setFileName(fnAll_);
 292         plotm->setTitle("Angle");
 293         plotm->setXAxisIsTime();
 294         plotm->setYLabel("Angle (degrees)");
 295         // TODO: Add legends? (there can be a massive amount of columns)
 296         angles_.addModule(plotm);
 297     }
 298 }
 299
 300
 301 //! Helper method to process selections into an array of coordinates.
 302 static void
 303 copy_pos(const SelectionList &sel, int natoms, int g, int first, rvec x[])
 304 {
 305     for (int k = 0; k < natoms; ++k)
 306     {
 307         copy_rvec(sel[g].position(first + k).x(), x[k]);
 308     }
 309 }
 310
 311
 312 //! Helper method to calculate a vector from two or three positions.
 313 static void
 314 calc_vec(int natoms, rvec x[], t_pbc *pbc, rvec xout, rvec cout)
 315 {
 316     switch (natoms)
 317     {
 318         case 2:
 319             if (pbc)
 320             {
 321                 pbc_dx(pbc, x[1], x[0], xout);
 322             }
 323             else
 324             {
 325                 rvec_sub(x[1], x[0], xout);
 326             }
 327             svmul(0.5, xout, cout);
 328             rvec_add(x[0], cout, cout);
 329             break;
 330         case 3:
 331         {
 332             rvec v1, v2;
 333             if (pbc)
 334             {
 335                 pbc_dx(pbc, x[1], x[0], v1);
 336                 pbc_dx(pbc, x[2], x[0], v2);
 337             }
 338             else
 339             {
 340                 rvec_sub(x[1], x[0], v1);
 341                 rvec_sub(x[2], x[0], v2);
 342             }
 343             cprod(v1, v2, xout);
 344             rvec_add(x[0], x[1], cout);
 345             rvec_add(cout, x[2], cout);
 346             svmul(1.0/3.0, cout, cout);
 347             break;
 348         }
 349         default:
 350             GMX_RELEASE_ASSERT(false, "Incorrectly initialized number of atoms");
 351     }
 352 }
 353
 354
 355 void
 356 Angle::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 357                     TrajectoryAnalysisModuleData *pdata)
 358 {
 359     AnalysisDataHandle       dh   = pdata->dataHandle(angles_);
 360     const SelectionList     &sel1 = pdata->parallelSelections(sel1_);
 361     const SelectionList     &sel2 = pdata->parallelSelections(sel2_);
 362
 363     checkSelections(sel1, sel2);
 364
 365     dh.startFrame(frnr, fr.time);
 366
 367     for (size_t g = 0; g < sel1_.size(); ++g)
 368     {
 369         rvec  v1, v2;
 370         rvec  c1, c2;
 371         switch (g2type_[0])
 372         {
 373             case 'z':
 374                 clear_rvec(v2);
 375                 v2[ZZ] = 1.0;
 376                 clear_rvec(c2);
 377                 break;
 378             case 's':
 379                 copy_rvec(sel2_[g].position(0).x(), c2);
 380                 break;
 381         }
 382         for (int i = 0, j = 0, n = 0;
 383              i < sel1[g].posCount();
 384              i += natoms1_, j += natoms2_, ++n)
 385         {
 386             rvec x[4];
 387             real angle;
 388             copy_pos(sel1, natoms1_, g, i, x);
 389             switch (g1type_[0])
 390             {
 391                 case 'a':
 392                     if (pbc)
 393                     {
 394                         pbc_dx(pbc, x[0], x[1], v1);
 395                         pbc_dx(pbc, x[2], x[1], v2);
 396                     }
 397                     else
 398                     {
 399                         rvec_sub(x[0], x[1], v1);
 400                         rvec_sub(x[2], x[1], v2);
 401                     }
 402                     angle = gmx_angle(v1, v2);
 403                     break;
 404                 case 'd':
 405                 {
 406                     rvec dx[3];
 407                     if (pbc)
 408                     {
 409                         pbc_dx(pbc, x[0], x[1], dx[0]);
 410                         pbc_dx(pbc, x[2], x[1], dx[1]);
 411                         pbc_dx(pbc, x[2], x[3], dx[2]);
 412                     }
 413                     else
 414                     {
 415                         rvec_sub(x[0], x[1], dx[0]);
 416                         rvec_sub(x[2], x[1], dx[1]);
 417                         rvec_sub(x[2], x[3], dx[2]);
 418                     }
 419                     cprod(dx[0], dx[1], v1);
 420                     cprod(dx[1], dx[2], v2);
 421                     angle = gmx_angle(v1, v2);
 422                     real ipr = iprod(dx[0], v2);
 423                     if (ipr < 0)
 424                     {
 425                         angle = -angle;
 426                     }
 427                     break;
 428                 }
 429                 case 'v':
 430                 case 'p':
 431                     calc_vec(natoms1_, x, pbc, v1, c1);
 432                     switch (g2type_[0])
 433                     {
 434                         case 'v':
 435                         case 'p':
 436                             copy_pos(sel2, natoms2_, 0, j, x);
 437                             calc_vec(natoms2_, x, pbc, v2, c2);
 438                             break;
 439                         case 't':
 440                             // FIXME: This is not parallelizable.
 441                             if (frnr == 0)
 442                             {
 443                                 copy_rvec(v1, vt0_[g][n]);
 444                             }
 445                             copy_rvec(vt0_[g][n], v2);
 446                             break;
 447                         case 'z':
 448                             c1[XX] = c1[YY] = 0.0;
 449                             break;
 450                         case 's':
 451                             if (pbc)
 452                             {
 453                                 pbc_dx(pbc, c1, c2, v2);
 454                             }
 455                             else
 456                             {
 457                                 rvec_sub(c1, c2, v2);
 458                             }
 459                             break;
 460                         default:
 461                             GMX_THROW(InternalError("invalid -g2 value"));
 462                     }
 463                     angle = gmx_angle(v1, v2);
 464                     break;
 465                 default:
 466                     GMX_THROW(InternalError("invalid -g1 value"));
 467             }
 468             /* TODO: Should we also calculate distances like g_sgangle?
 469              * Could be better to leave that for a separate tool.
 470                real dist = 0.0;
 471                if (bDumpDist_)
 472                {
 473                 if (pbc)
 474                 {
 475                     rvec dx;
 476                     pbc_dx(pbc, c2, c1, dx);
 477                     dist = norm(dx);
 478                 }
 479                 else
 480                 {
 481                     dist = sqrt(distance2(c1, c2));
 482                 }
 483                }
 484              */
 485             dh.setPoint(n, angle * RAD2DEG);
 486         }
 487     }
 488     dh.finishFrame();
 489 }
 490
 491
 492 void
 493 Angle::finishAnalysis(int /*nframes*/)
 494 {
 495 }
 496
 497
 498 void
 499 Angle::writeOutput()
 500 {
 501 }
 502
 503 } // namespace modules
 504
 505 } // namespace gmxana