src/gromacs/trajectoryanalysis/modules/select.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2009,2010,2011,2012,2013 by the GROMACS development team.
   5  * Copyright (c) 2014,2015,2016,2017,2018 by the GROMACS development team.
   6  * Copyright (c) 2019,2020, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 /*! \internal \file
  38  * \brief
  39  * Implements gmx::analysismodules::Select.
  40  *
  41  * \author Teemu Murtola <teemu.murtola@gmail.com>
  42  * \ingroup module_trajectoryanalysis
  43  */
  44 #include "gmxpre.h"
  45
  46 #include "select.h"
  47
  48 #include <cstdio>
  49 #include <cstring>
  50
  51 #include <algorithm>
  52 #include <set>
  53 #include <string>
  54 #include <vector>
  55
  56 #include "gromacs/analysisdata/analysisdata.h"
  57 #include "gromacs/analysisdata/dataframe.h"
  58 #include "gromacs/analysisdata/datamodule.h"
  59 #include "gromacs/analysisdata/modules/average.h"
  60 #include "gromacs/analysisdata/modules/lifetime.h"
  61 #include "gromacs/analysisdata/modules/plot.h"
  62 #include "gromacs/fileio/gmxfio.h"
  63 #include "gromacs/fileio/trxio.h"
  64 #include "gromacs/options/basicoptions.h"
  65 #include "gromacs/options/filenameoption.h"
  66 #include "gromacs/options/ioptionscontainer.h"
  67 #include "gromacs/selection/selection.h"
  68 #include "gromacs/selection/selectionoption.h"
  69 #include "gromacs/topology/topology.h"
  70 #include "gromacs/trajectory/trajectoryframe.h"
  71 #include "gromacs/trajectoryanalysis/analysissettings.h"
  72 #include "gromacs/trajectoryanalysis/topologyinformation.h"
  73 #include "gromacs/utility/arrayref.h"
  74 #include "gromacs/utility/exceptions.h"
  75 #include "gromacs/utility/gmxassert.h"
  76 #include "gromacs/utility/smalloc.h"
  77 #include "gromacs/utility/stringutil.h"
  78 #include "gromacs/utility/unique_cptr.h"
  79
  80 namespace gmx
  81 {
  82
  83 namespace analysismodules
  84 {
  85
  86 namespace
  87 {
  88
  89 /*! \internal \brief
  90  * Data module for writing index files.
  91  *
  92  * \ingroup module_analysisdata
  93  */
  94 class IndexFileWriterModule : public AnalysisDataModuleSerial
  95 {
  96 public:
  97     IndexFileWriterModule();
  98     ~IndexFileWriterModule() override;
  99
 100     //! Sets the file name to write the index file to.
 101     void setFileName(const std::string& fnm);
 102     /*! \brief
 103      * Adds information about a group to be printed.
 104      *
 105      * Must be called for each group present in the input data.
 106      */
 107     void addGroup(const std::string& name, bool bDynamic);
 108
 109     int flags() const override;
 110
 111     void dataStarted(AbstractAnalysisData* data) override;
 112     void frameStarted(const AnalysisDataFrameHeader& header) override;
 113     void pointsAdded(const AnalysisDataPointSetRef& points) override;
 114     void frameFinished(const AnalysisDataFrameHeader& header) override;
 115     void dataFinished() override;
 116
 117 private:
 118     void closeFile();
 119
 120     struct GroupInfo
 121     {
 122         GroupInfo(const std::string& name, bool bDynamic) : name(name), bDynamic(bDynamic) {}
 123
 124         std::string name;
 125         bool        bDynamic;
 126     };
 127
 128     std::string            fnm_;
 129     std::vector<GroupInfo> groups_;
 130     FILE*                  fp_;
 131     int                    currentGroup_;
 132     int                    currentSize_;
 133     bool                   bAnyWritten_;
 134 };
 135
 136 /********************************************************************
 137  * IndexFileWriterModule
 138  */
 139
 140 IndexFileWriterModule::IndexFileWriterModule() :
 141     fp_(nullptr),
 142     currentGroup_(-1),
 143     currentSize_(0),
 144     bAnyWritten_(false)
 145 {
 146 }
 147
 148
 149 IndexFileWriterModule::~IndexFileWriterModule()
 150 {
 151     closeFile();
 152 }
 153
 154
 155 void IndexFileWriterModule::closeFile()
 156 {
 157     if (fp_ != nullptr)
 158     {
 159         gmx_fio_fclose(fp_);
 160         fp_ = nullptr;
 161     }
 162 }
 163
 164
 165 void IndexFileWriterModule::setFileName(const std::string& fnm)
 166 {
 167     fnm_ = fnm;
 168 }
 169
 170
 171 void IndexFileWriterModule::addGroup(const std::string& name, bool bDynamic)
 172 {
 173     std::string newName(name);
 174     std::replace(newName.begin(), newName.end(), ' ', '_');
 175     groups_.emplace_back(newName, bDynamic);
 176 }
 177
 178
 179 int IndexFileWriterModule::flags() const
 180 {
 181     return efAllowMulticolumn | efAllowMultipoint;
 182 }
 183
 184
 185 void IndexFileWriterModule::dataStarted(AbstractAnalysisData* /*data*/)
 186 {
 187     if (!fnm_.empty())
 188     {
 189         fp_ = gmx_fio_fopen(fnm_.c_str(), "w");
 190     }
 191 }
 192
 193
 194 void IndexFileWriterModule::frameStarted(const AnalysisDataFrameHeader& /*header*/)
 195 {
 196     bAnyWritten_  = false;
 197     currentGroup_ = -1;
 198 }
 199
 200
 201 void IndexFileWriterModule::pointsAdded(const AnalysisDataPointSetRef& points)
 202 {
 203     if (fp_ == nullptr)
 204     {
 205         return;
 206     }
 207     bool bFirstFrame = (points.frameIndex() == 0);
 208     if (points.firstColumn() == 0)
 209     {
 210         ++currentGroup_;
 211         GMX_RELEASE_ASSERT(currentGroup_ < ssize(groups_), "Too few groups initialized");
 212         if (bFirstFrame || groups_[currentGroup_].bDynamic)
 213         {
 214             if (!bFirstFrame || currentGroup_ > 0)
 215             {
 216                 std::fprintf(fp_, "\n\n");
 217             }
 218             std::string name = groups_[currentGroup_].name;
 219             if (groups_[currentGroup_].bDynamic)
 220             {
 221                 name += formatString("_f%d_t%.3f", points.frameIndex(), points.x());
 222             }
 223             std::fprintf(fp_, "[ %s ]", name.c_str());
 224             bAnyWritten_ = true;
 225             currentSize_ = 0;
 226         }
 227     }
 228     else
 229     {
 230         if (bFirstFrame || groups_[currentGroup_].bDynamic)
 231         {
 232             if (currentSize_ % 15 == 0)
 233             {
 234                 std::fprintf(fp_, "\n");
 235             }
 236             std::fprintf(fp_, "%4d ", static_cast<int>(points.y(0)));
 237             ++currentSize_;
 238         }
 239     }
 240 }
 241
 242
 243 void IndexFileWriterModule::frameFinished(const AnalysisDataFrameHeader& /*header*/) {}
 244
 245
 246 void IndexFileWriterModule::dataFinished()
 247 {
 248     if (fp_ != nullptr)
 249     {
 250         std::fprintf(fp_, "\n");
 251     }
 252     closeFile();
 253 }
 254
 255 /********************************************************************
 256  * Select
 257  */
 258
 259 //! How to identify residues in output files.
 260 enum class ResidueNumbering : int
 261 {
 262     ByNumber,
 263     ByIndex,
 264     Count
 265 };
 266 //! Which atoms to write out to PDB files.
 267 enum class PdbAtomsSelection : int
 268 {
 269     All,
 270     MaxSelection,
 271     Selected,
 272     Count
 273 };
 274 //! String values corresponding to ResidueNumbering.
 275 const EnumerationArray<ResidueNumbering, const char*> c_residueNumberingTypeNames = { { "number",
 276                                                                                         "index" } };
 277 //! String values corresponding to PdbAtomsSelection.
 278 const EnumerationArray<PdbAtomsSelection, const char*> c_pdbAtomsTypeNames = { { "all", "maxsel",
 279                                                                                  "selected" } };
 280
 281 class Select : public TrajectoryAnalysisModule
 282 {
 283 public:
 284     Select();
 285
 286     void initOptions(IOptionsContainer* options, TrajectoryAnalysisSettings* settings) override;
 287     void optionsFinished(TrajectoryAnalysisSettings* settings) override;
 288     void initAnalysis(const TrajectoryAnalysisSettings& settings, const TopologyInformation& top) override;
 289
 290     void analyzeFrame(int frnr, const t_trxframe& fr, t_pbc* pbc, TrajectoryAnalysisModuleData* pdata) override;
 291
 292     void finishAnalysis(int nframes) override;
 293     void writeOutput() override;
 294
 295 private:
 296     SelectionList sel_;
 297
 298     std::string       fnSize_;
 299     std::string       fnFrac_;
 300     std::string       fnIndex_;
 301     std::string       fnNdx_;
 302     std::string       fnMask_;
 303     std::string       fnOccupancy_;
 304     std::string       fnPDB_;
 305     std::string       fnLifetime_;
 306     bool              bTotNorm_;
 307     bool              bFracNorm_;
 308     bool              bResInd_;
 309     bool              bCumulativeLifetimes_;
 310     ResidueNumbering  resNumberType_;
 311     PdbAtomsSelection pdbAtoms_;
 312
 313     //! The input topology.
 314     const TopologyInformation*        top_;
 315     std::vector<int>                  totsize_;
 316     AnalysisData                      sdata_;
 317     AnalysisData                      cdata_;
 318     AnalysisData                      idata_;
 319     AnalysisData                      mdata_;
 320     AnalysisDataAverageModulePointer  occupancyModule_;
 321     AnalysisDataLifetimeModulePointer lifetimeModule_;
 322 };
 323
 324 Select::Select() :
 325     bTotNorm_(false),
 326     bFracNorm_(false),
 327     bResInd_(false),
 328     bCumulativeLifetimes_(true),
 329     resNumberType_(ResidueNumbering::ByNumber),
 330     pdbAtoms_(PdbAtomsSelection::All),
 331     top_(nullptr),
 332     occupancyModule_(new AnalysisDataAverageModule()),
 333     lifetimeModule_(new AnalysisDataLifetimeModule())
 334 {
 335     mdata_.addModule(occupancyModule_);
 336     mdata_.addModule(lifetimeModule_);
 337
 338     registerAnalysisDataset(&sdata_, "size");
 339     registerAnalysisDataset(&cdata_, "cfrac");
 340     idata_.setColumnCount(0, 2);
 341     idata_.setMultipoint(true);
 342     registerAnalysisDataset(&idata_, "index");
 343     registerAnalysisDataset(&mdata_, "mask");
 344     occupancyModule_->setXAxis(1.0, 1.0);
 345     registerBasicDataset(occupancyModule_.get(), "occupancy");
 346     registerBasicDataset(lifetimeModule_.get(), "lifetime");
 347 }
 348
 349
 350 void Select::initOptions(IOptionsContainer* options, TrajectoryAnalysisSettings* settings)
 351 {
 352     static const char* const desc[] = {
 353         "[THISMODULE] writes out basic data about dynamic selections.",
 354         "It can be used for some simple analyses, or the output can",
 355         "be combined with output from other programs and/or external",
 356         "analysis programs to calculate more complex things.",
 357         "For detailed help on the selection syntax, please use",
 358         "[TT]gmx help selections[tt].",
 359         "",
 360         "Any combination of the output options is possible, but note",
 361         "that [TT]-om[tt] only operates on the first selection.",
 362         "Also note that if you provide no output options, no output is",
 363         "produced.",
 364         "",
 365         "With [TT]-os[tt], calculates the number of positions in each",
 366         "selection for each frame. With [TT]-norm[tt], the output is",
 367         "between 0 and 1 and describes the fraction from the maximum",
 368         "number of positions (e.g., for selection 'resname RA and x < 5'",
 369         "the maximum number of positions is the number of atoms in",
 370         "RA residues). With [TT]-cfnorm[tt], the output is divided",
 371         "by the fraction covered by the selection.",
 372         "[TT]-norm[tt] and [TT]-cfnorm[tt] can be specified independently",
 373         "of one another.",
 374         "",
 375         "With [TT]-oc[tt], the fraction covered by each selection is",
 376         "written out as a function of time.",
 377         "",
 378         "With [TT]-oi[tt], the selected atoms/residues/molecules are",
 379         "written out as a function of time. In the output, the first",
 380         "column contains the frame time, the second contains the number",
 381         "of positions, followed by the atom/residue/molecule numbers.",
 382         "If more than one selection is specified, the size of the second",
 383         "group immediately follows the last number of the first group",
 384         "and so on.",
 385         "",
 386         "With [TT]-on[tt], the selected atoms are written as a index file",
 387         "compatible with [TT]make_ndx[tt] and the analyzing tools. Each selection",
 388         "is written as a selection group and for dynamic selections a",
 389         "group is written for each frame.",
 390         "",
 391         "For residue numbers, the output of [TT]-oi[tt] can be controlled",
 392         "with [TT]-resnr[tt]: [TT]number[tt] (default) prints the residue",
 393         "numbers as they appear in the input file, while [TT]index[tt] prints",
 394         "unique numbers assigned to the residues in the order they appear",
 395         "in the input file, starting with 1. The former is more intuitive,",
 396         "but if the input contains multiple residues with the same number,",
 397         "the output can be less useful.",
 398         "",
 399         "With [TT]-om[tt], a mask is printed for the first selection",
 400         "as a function of time. Each line in the output corresponds to",
 401         "one frame, and contains either 0/1 for each atom/residue/molecule",
 402         "possibly selected. 1 stands for the atom/residue/molecule being",
 403         "selected for the current frame, 0 for not selected.",
 404         "",
 405         "With [TT]-of[tt], the occupancy fraction of each position (i.e.,",
 406         "the fraction of frames where the position is selected) is",
 407         "printed.",
 408         "",
 409         "With [TT]-ofpdb[tt], a PDB file is written out where the occupancy",
 410         "column is filled with the occupancy fraction of each atom in the",
 411         "selection. The coordinates in the PDB file will be those from the",
 412         "input topology. [TT]-pdbatoms[tt] can be used to control which atoms",
 413         "appear in the output PDB file: with [TT]all[tt] all atoms are",
 414         "present, with [TT]maxsel[tt] all atoms possibly selected by the",
 415         "selection are present, and with [TT]selected[tt] only atoms that are",
 416         "selected at least in one frame are present.",
 417         "",
 418         "With [TT]-olt[tt], a histogram is produced that shows the number of",
 419         "selected positions as a function of the time the position was",
 420         "continuously selected. [TT]-cumlt[tt] can be used to control whether",
 421         "subintervals of longer intervals are included in the histogram.[PAR]",
 422         "[TT]-om[tt], [TT]-of[tt], and [TT]-olt[tt] only make sense with",
 423         "dynamic selections.",
 424         "",
 425         "To plot coordinates for selections, use [gmx-trajectory]."
 426     };
 427
 428     settings->setHelpText(desc);
 429
 430     options->addOption(FileNameOption("os")
 431                                .filetype(eftPlot)
 432                                .outputFile()
 433                                .store(&fnSize_)
 434                                .defaultBasename("size")
 435                                .description("Number of positions in each selection"));
 436     options->addOption(FileNameOption("oc")
 437                                .filetype(eftPlot)
 438                                .outputFile()
 439                                .store(&fnFrac_)
 440                                .defaultBasename("cfrac")
 441                                .description("Covered fraction for each selection"));
 442     options->addOption(FileNameOption("oi")
 443                                .filetype(eftGenericData)
 444                                .outputFile()
 445                                .store(&fnIndex_)
 446                                .defaultBasename("index")
 447                                .description("Indices selected by each selection"));
 448     options->addOption(FileNameOption("on")
 449                                .filetype(eftIndex)
 450                                .outputFile()
 451                                .store(&fnNdx_)
 452                                .defaultBasename("index")
 453                                .description("Index file from the selection"));
 454     options->addOption(FileNameOption("om")
 455                                .filetype(eftPlot)
 456                                .outputFile()
 457                                .store(&fnMask_)
 458                                .defaultBasename("mask")
 459                                .description("Mask for selected positions"));
 460     options->addOption(FileNameOption("of")
 461                                .filetype(eftPlot)
 462                                .outputFile()
 463                                .store(&fnOccupancy_)
 464                                .defaultBasename("occupancy")
 465                                .description("Occupied fraction for selected positions"));
 466     options->addOption(
 467             FileNameOption("ofpdb")
 468                     .filetype(eftPDB)
 469                     .outputFile()
 470                     .store(&fnPDB_)
 471                     .defaultBasename("occupancy")
 472                     .description("PDB file with occupied fraction for selected positions"));
 473     options->addOption(FileNameOption("olt")
 474                                .filetype(eftPlot)
 475                                .outputFile()
 476                                .store(&fnLifetime_)
 477                                .defaultBasename("lifetime")
 478                                .description("Lifetime histogram"));
 479
 480     options->addOption(SelectionOption("select").storeVector(&sel_).required().multiValue().description(
 481             "Selections to analyze"));
 482
 483     options->addOption(
 484             BooleanOption("norm").store(&bTotNorm_).description("Normalize by total number of positions with -os"));
 485     options->addOption(
 486             BooleanOption("cfnorm").store(&bFracNorm_).description("Normalize by covered fraction with -os"));
 487     options->addOption(EnumOption<ResidueNumbering>("resnr")
 488                                .store(&resNumberType_)
 489                                .enumValue(c_residueNumberingTypeNames)
 490                                .description("Residue number output type with -oi and -on"));
 491     options->addOption(EnumOption<PdbAtomsSelection>("pdbatoms")
 492                                .store(&pdbAtoms_)
 493                                .enumValue(c_pdbAtomsTypeNames)
 494                                .description("Atoms to write with -ofpdb"));
 495     options->addOption(BooleanOption("cumlt")
 496                                .store(&bCumulativeLifetimes_)
 497                                .description("Cumulate subintervals of longer intervals in -olt"));
 498 }
 499
 500 void Select::optionsFinished(TrajectoryAnalysisSettings* settings)
 501 {
 502     if (!fnPDB_.empty())
 503     {
 504         settings->setFlag(TrajectoryAnalysisSettings::efRequireTop);
 505         settings->setFlag(TrajectoryAnalysisSettings::efUseTopX);
 506     }
 507 }
 508
 509 void Select::initAnalysis(const TrajectoryAnalysisSettings& settings, const TopologyInformation& top)
 510 {
 511     bResInd_ = (resNumberType_ == ResidueNumbering::ByIndex);
 512
 513     for (SelectionList::iterator i = sel_.begin(); i != sel_.end(); ++i)
 514     {
 515         i->initCoveredFraction(CFRAC_SOLIDANGLE);
 516     }
 517
 518     // TODO: For large systems, a float may not have enough precision
 519     sdata_.setColumnCount(0, sel_.size());
 520     totsize_.reserve(sel_.size());
 521     for (size_t g = 0; g < sel_.size(); ++g)
 522     {
 523         totsize_.push_back(sel_[g].posCount());
 524     }
 525     if (!fnSize_.empty())
 526     {
 527         AnalysisDataPlotModulePointer plot(new AnalysisDataPlotModule(settings.plotSettings()));
 528         plot->setFileName(fnSize_);
 529         plot->setTitle("Selection size");
 530         plot->setXAxisIsTime();
 531         plot->setYLabel("Number");
 532         for (size_t g = 0; g < sel_.size(); ++g)
 533         {
 534             plot->appendLegend(sel_[g].name());
 535         }
 536         sdata_.addModule(plot);
 537     }
 538
 539     cdata_.setColumnCount(0, sel_.size());
 540     if (!fnFrac_.empty())
 541     {
 542         AnalysisDataPlotModulePointer plot(new AnalysisDataPlotModule(settings.plotSettings()));
 543         plot->setFileName(fnFrac_);
 544         plot->setTitle("Covered fraction");
 545         plot->setXAxisIsTime();
 546         plot->setYLabel("Fraction");
 547         plot->setYFormat(6, 4);
 548         for (size_t g = 0; g < sel_.size(); ++g)
 549         {
 550             plot->appendLegend(sel_[g].name());
 551         }
 552         cdata_.addModule(plot);
 553     }
 554
 555     // TODO: For large systems, a float may not have enough precision
 556     if (!fnIndex_.empty())
 557     {
 558         AnalysisDataPlotModulePointer plot(new AnalysisDataPlotModule(settings.plotSettings()));
 559         plot->setFileName(fnIndex_);
 560         plot->setPlainOutput(true);
 561         plot->setYFormat(4, 0);
 562         idata_.addModule(plot);
 563     }
 564     if (!fnNdx_.empty())
 565     {
 566         std::shared_ptr<IndexFileWriterModule> writer(new IndexFileWriterModule());
 567         writer->setFileName(fnNdx_);
 568         for (size_t g = 0; g < sel_.size(); ++g)
 569         {
 570             writer->addGroup(sel_[g].name(), sel_[g].isDynamic());
 571         }
 572         idata_.addModule(writer);
 573     }
 574
 575     mdata_.setDataSetCount(sel_.size());
 576     for (size_t g = 0; g < sel_.size(); ++g)
 577     {
 578         mdata_.setColumnCount(g, sel_[g].posCount());
 579     }
 580     lifetimeModule_->setCumulative(bCumulativeLifetimes_);
 581     if (!fnMask_.empty())
 582     {
 583         AnalysisDataPlotModulePointer plot(new AnalysisDataPlotModule(settings.plotSettings()));
 584         plot->setFileName(fnMask_);
 585         plot->setTitle("Selection mask");
 586         plot->setXAxisIsTime();
 587         plot->setYLabel("Occupancy");
 588         plot->setYFormat(1, 0);
 589         // TODO: Add legend? (there can be massive amount of columns)
 590         mdata_.addModule(plot);
 591     }
 592     if (!fnOccupancy_.empty())
 593     {
 594         AnalysisDataPlotModulePointer plot(new AnalysisDataPlotModule(settings.plotSettings()));
 595         plot->setFileName(fnOccupancy_);
 596         plot->setTitle("Fraction of time selection matches");
 597         plot->setXLabel("Selected position");
 598         plot->setYLabel("Occupied fraction");
 599         for (size_t g = 0; g < sel_.size(); ++g)
 600         {
 601             plot->appendLegend(sel_[g].name());
 602         }
 603         occupancyModule_->addModule(plot);
 604     }
 605     if (!fnLifetime_.empty())
 606     {
 607         AnalysisDataPlotModulePointer plot(new AnalysisDataPlotModule(settings.plotSettings()));
 608         plot->setFileName(fnLifetime_);
 609         plot->setTitle("Lifetime histogram");
 610         plot->setXAxisIsTime();
 611         plot->setYLabel("Number of occurrences");
 612         for (size_t g = 0; g < sel_.size(); ++g)
 613         {
 614             plot->appendLegend(sel_[g].name());
 615         }
 616         lifetimeModule_->addModule(plot);
 617     }
 618
 619     top_ = &top;
 620 }
 621
 622
 623 void Select::analyzeFrame(int frnr, const t_trxframe& fr, t_pbc* /* pbc */, TrajectoryAnalysisModuleData* pdata)
 624 {
 625     AnalysisDataHandle   sdh = pdata->dataHandle(sdata_);
 626     AnalysisDataHandle   cdh = pdata->dataHandle(cdata_);
 627     AnalysisDataHandle   idh = pdata->dataHandle(idata_);
 628     AnalysisDataHandle   mdh = pdata->dataHandle(mdata_);
 629     const SelectionList& sel = pdata->parallelSelections(sel_);
 630
 631     sdh.startFrame(frnr, fr.time);
 632     for (size_t g = 0; g < sel.size(); ++g)
 633     {
 634         real normfac = bFracNorm_ ? 1.0 / sel[g].coveredFraction() : 1.0;
 635         if (bTotNorm_)
 636         {
 637             normfac /= totsize_[g];
 638         }
 639         sdh.setPoint(g, sel[g].posCount() * normfac);
 640     }
 641     sdh.finishFrame();
 642
 643     cdh.startFrame(frnr, fr.time);
 644     for (size_t g = 0; g < sel.size(); ++g)
 645     {
 646         cdh.setPoint(g, sel[g].coveredFraction());
 647     }
 648     cdh.finishFrame();
 649
 650     idh.startFrame(frnr, fr.time);
 651     for (size_t g = 0; g < sel.size(); ++g)
 652     {
 653         idh.setPoint(0, sel[g].posCount());
 654         idh.finishPointSet();
 655         for (int i = 0; i < sel[g].posCount(); ++i)
 656         {
 657             const SelectionPosition& p = sel[g].position(i);
 658             if (sel[g].type() == INDEX_RES && !bResInd_)
 659             {
 660                 idh.setPoint(1, top_->atoms()->resinfo[p.mappedId()].nr);
 661             }
 662             else
 663             {
 664                 idh.setPoint(1, p.mappedId() + 1);
 665             }
 666             idh.finishPointSet();
 667         }
 668     }
 669     idh.finishFrame();
 670
 671     mdh.startFrame(frnr, fr.time);
 672     for (size_t g = 0; g < sel.size(); ++g)
 673     {
 674         mdh.selectDataSet(g);
 675         for (int i = 0; i < totsize_[g]; ++i)
 676         {
 677             mdh.setPoint(i, 0);
 678         }
 679         for (int i = 0; i < sel[g].posCount(); ++i)
 680         {
 681             mdh.setPoint(sel[g].position(i).refId(), 1);
 682         }
 683     }
 684     mdh.finishFrame();
 685 }
 686
 687
 688 void Select::finishAnalysis(int /*nframes*/) {}
 689
 690
 691 void Select::writeOutput()
 692 {
 693     if (!fnPDB_.empty())
 694     {
 695         GMX_RELEASE_ASSERT(top_->hasTopology(),
 696                            "Topology should have been loaded or an error given earlier");
 697         auto atoms = top_->copyAtoms();
 698         if (!atoms->havePdbInfo)
 699         {
 700             snew(atoms->pdbinfo, atoms->nr);
 701             atoms->havePdbInfo = TRUE;
 702         }
 703         for (int i = 0; i < atoms->nr; ++i)
 704         {
 705             atoms->pdbinfo[i].occup = 0.0;
 706         }
 707         for (size_t g = 0; g < sel_.size(); ++g)
 708         {
 709             for (int i = 0; i < sel_[g].posCount(); ++i)
 710             {
 711                 ArrayRef<const int>                 atomIndices = sel_[g].position(i).atomIndices();
 712                 ArrayRef<const int>::const_iterator ai;
 713                 for (ai = atomIndices.begin(); ai != atomIndices.end(); ++ai)
 714                 {
 715                     atoms->pdbinfo[*ai].occup += occupancyModule_->average(g, i);
 716                 }
 717             }
 718         }
 719
 720         std::vector<RVec> x = copyOf(top_->x());
 721         t_trxframe        fr;
 722         clear_trxframe(&fr, TRUE);
 723         fr.bAtoms = TRUE;
 724         fr.atoms  = atoms.get();
 725         fr.bX     = TRUE;
 726         fr.bBox   = TRUE;
 727         fr.x      = as_rvec_array(x.data());
 728         top_->getBox(fr.box);
 729
 730         switch (pdbAtoms_)
 731         {
 732             case PdbAtomsSelection::All:
 733             {
 734                 t_trxstatus* status = open_trx(fnPDB_.c_str(), "w");
 735                 write_trxframe(status, &fr, nullptr);
 736                 close_trx(status);
 737                 break;
 738             }
 739             case PdbAtomsSelection::MaxSelection:
 740             {
 741                 std::set<int> atomIndicesSet;
 742                 for (size_t g = 0; g < sel_.size(); ++g)
 743                 {
 744                     ArrayRef<const int> atomIndices = sel_[g].atomIndices();
 745                     atomIndicesSet.insert(atomIndices.begin(), atomIndices.end());
 746                 }
 747                 std::vector<int> allAtomIndices(atomIndicesSet.begin(), atomIndicesSet.end());
 748                 t_trxstatus*     status = open_trx(fnPDB_.c_str(), "w");
 749                 write_trxframe_indexed(status, &fr, allAtomIndices.size(), allAtomIndices.data(), nullptr);
 750                 close_trx(status);
 751                 break;
 752             }
 753             case PdbAtomsSelection::Selected:
 754             {
 755                 std::vector<int> indices;
 756                 for (int i = 0; i < atoms->nr; ++i)
 757                 {
 758                     if (atoms->pdbinfo[i].occup > 0.0)
 759                     {
 760                         indices.push_back(i);
 761                     }
 762                 }
 763                 t_trxstatus* status = open_trx(fnPDB_.c_str(), "w");
 764                 write_trxframe_indexed(status, &fr, indices.size(), indices.data(), nullptr);
 765                 close_trx(status);
 766                 break;
 767             }
 768             case PdbAtomsSelection::Count:
 769                 GMX_RELEASE_ASSERT(false,
 770                                    "Mismatch between -pdbatoms enum values and implementation");
 771         }
 772     }
 773 }
 774
 775 } // namespace
 776
 777 const char SelectInfo::name[]             = "select";
 778 const char SelectInfo::shortDescription[] = "Print general information about selections";
 779
 780 TrajectoryAnalysisModulePointer SelectInfo::create()
 781 {
 782     return TrajectoryAnalysisModulePointer(new Select);
 783 }
 784
 785 } // namespace analysismodules
 786
 787 } // namespace gmx