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