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