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,2015, 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/ioptionscontainer.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(IOptionsContainer          *options,
 265                                  TrajectoryAnalysisSettings *settings);
 266         virtual void optionsFinished(TrajectoryAnalysisSettings *settings);
 267         virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
 268                                   const TopologyInformation        &top);
 269
 270         virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 271                                   TrajectoryAnalysisModuleData *pdata);
 272
 273         virtual void finishAnalysis(int nframes);
 274         virtual void writeOutput();
 275
 276     private:
 277         SelectionList                       sel_;
 278         SelectionOptionInfo                *selOpt_;
 279
 280         std::string                         fnSize_;
 281         std::string                         fnFrac_;
 282         std::string                         fnIndex_;
 283         std::string                         fnNdx_;
 284         std::string                         fnMask_;
 285         std::string                         fnOccupancy_;
 286         std::string                         fnPDB_;
 287         std::string                         fnLifetime_;
 288         bool                                bTotNorm_;
 289         bool                                bFracNorm_;
 290         bool                                bResInd_;
 291         bool                                bCumulativeLifetimes_;
 292         std::string                         resNumberType_;
 293         std::string                         pdbAtoms_;
 294
 295         const TopologyInformation          *top_;
 296         std::vector<int>                    totsize_;
 297         AnalysisData                        sdata_;
 298         AnalysisData                        cdata_;
 299         AnalysisData                        idata_;
 300         AnalysisData                        mdata_;
 301         AnalysisDataAverageModulePointer    occupancyModule_;
 302         AnalysisDataLifetimeModulePointer   lifetimeModule_;
 303 };
 304
 305 Select::Select()
 306     : TrajectoryAnalysisModule(SelectInfo::name, SelectInfo::shortDescription),
 307       selOpt_(NULL),
 308       bTotNorm_(false), bFracNorm_(false), bResInd_(false),
 309       bCumulativeLifetimes_(true), top_(NULL),
 310       occupancyModule_(new AnalysisDataAverageModule()),
 311       lifetimeModule_(new AnalysisDataLifetimeModule())
 312 {
 313     mdata_.addModule(occupancyModule_);
 314     mdata_.addModule(lifetimeModule_);
 315
 316     registerAnalysisDataset(&sdata_, "size");
 317     registerAnalysisDataset(&cdata_, "cfrac");
 318     idata_.setColumnCount(0, 2);
 319     idata_.setMultipoint(true);
 320     registerAnalysisDataset(&idata_, "index");
 321     registerAnalysisDataset(&mdata_, "mask");
 322     occupancyModule_->setXAxis(1.0, 1.0);
 323     registerBasicDataset(occupancyModule_.get(), "occupancy");
 324     registerBasicDataset(lifetimeModule_.get(), "lifetime");
 325 }
 326
 327
 328 void
 329 Select::initOptions(IOptionsContainer *options, TrajectoryAnalysisSettings *settings)
 330 {
 331     static const char *const desc[] = {
 332         "[THISMODULE] writes out basic data about dynamic selections.",
 333         "It can be used for some simple analyses, or the output can",
 334         "be combined with output from other programs and/or external",
 335         "analysis programs to calculate more complex things.",
 336         "Any combination of the output options is possible, but note",
 337         "that [TT]-om[tt] only operates on the first selection.",
 338         "Also note that if you provide no output options, no output is",
 339         "produced.[PAR]",
 340         "With [TT]-os[tt], calculates the number of positions in each",
 341         "selection for each frame. With [TT]-norm[tt], the output is",
 342         "between 0 and 1 and describes the fraction from the maximum",
 343         "number of positions (e.g., for selection 'resname RA and x < 5'",
 344         "the maximum number of positions is the number of atoms in",
 345         "RA residues). With [TT]-cfnorm[tt], the output is divided",
 346         "by the fraction covered by the selection.",
 347         "[TT]-norm[tt] and [TT]-cfnorm[tt] can be specified independently",
 348         "of one another.[PAR]",
 349         "With [TT]-oc[tt], the fraction covered by each selection is",
 350         "written out as a function of time.[PAR]",
 351         "With [TT]-oi[tt], the selected atoms/residues/molecules are",
 352         "written out as a function of time. In the output, the first",
 353         "column contains the frame time, the second contains the number",
 354         "of positions, followed by the atom/residue/molecule numbers.",
 355         "If more than one selection is specified, the size of the second",
 356         "group immediately follows the last number of the first group",
 357         "and so on.[PAR]",
 358         "With [TT]-on[tt], the selected atoms are written as a index file",
 359         "compatible with [TT]make_ndx[tt] and the analyzing tools. Each selection",
 360         "is written as a selection group and for dynamic selections a",
 361         "group is written for each frame.[PAR]",
 362         "For residue numbers, the output of [TT]-oi[tt] can be controlled",
 363         "with [TT]-resnr[tt]: [TT]number[tt] (default) prints the residue",
 364         "numbers as they appear in the input file, while [TT]index[tt] prints",
 365         "unique numbers assigned to the residues in the order they appear",
 366         "in the input file, starting with 1. The former is more intuitive,",
 367         "but if the input contains multiple residues with the same number,",
 368         "the output can be less useful.[PAR]",
 369         "With [TT]-om[tt], a mask is printed for the first selection",
 370         "as a function of time. Each line in the output corresponds to",
 371         "one frame, and contains either 0/1 for each atom/residue/molecule",
 372         "possibly selected. 1 stands for the atom/residue/molecule being",
 373         "selected for the current frame, 0 for not selected.[PAR]",
 374         "With [TT]-of[tt], the occupancy fraction of each position (i.e.,",
 375         "the fraction of frames where the position is selected) is",
 376         "printed.[PAR]",
 377         "With [TT]-ofpdb[tt], a PDB file is written out where the occupancy",
 378         "column is filled with the occupancy fraction of each atom in the",
 379         "selection. The coordinates in the PDB file will be those from the",
 380         "input topology. [TT]-pdbatoms[tt] can be used to control which atoms",
 381         "appear in the output PDB file: with [TT]all[tt] all atoms are",
 382         "present, with [TT]maxsel[tt] all atoms possibly selected by the",
 383         "selection are present, and with [TT]selected[tt] only atoms that are",
 384         "selected at least in one frame are present.[PAR]",
 385         "With [TT]-olt[tt], a histogram is produced that shows the number of",
 386         "selected positions as a function of the time the position was",
 387         "continuously selected. [TT]-cumlt[tt] can be used to control whether",
 388         "subintervals of longer intervals are included in the histogram.[PAR]",
 389         "[TT]-om[tt], [TT]-of[tt], and [TT]-olt[tt] only make sense with",
 390         "dynamic selections."
 391     };
 392
 393     settings->setHelpText(desc);
 394
 395     options->addOption(FileNameOption("os").filetype(eftPlot).outputFile()
 396                            .store(&fnSize_).defaultBasename("size")
 397                            .description("Number of positions in each selection"));
 398     options->addOption(FileNameOption("oc").filetype(eftPlot).outputFile()
 399                            .store(&fnFrac_).defaultBasename("cfrac")
 400                            .description("Covered fraction for each selection"));
 401     options->addOption(FileNameOption("oi").filetype(eftGenericData).outputFile()
 402                            .store(&fnIndex_).defaultBasename("index")
 403                            .description("Indices selected by each selection"));
 404     options->addOption(FileNameOption("on").filetype(eftIndex).outputFile()
 405                            .store(&fnNdx_).defaultBasename("index")
 406                            .description("Index file from the selection"));
 407     options->addOption(FileNameOption("om").filetype(eftPlot).outputFile()
 408                            .store(&fnMask_).defaultBasename("mask")
 409                            .description("Mask for selected positions"));
 410     options->addOption(FileNameOption("of").filetype(eftPlot).outputFile()
 411                            .store(&fnOccupancy_).defaultBasename("occupancy")
 412                            .description("Occupied fraction for selected positions"));
 413     options->addOption(FileNameOption("ofpdb").filetype(eftPDB).outputFile()
 414                            .store(&fnPDB_).defaultBasename("occupancy")
 415                            .description("PDB file with occupied fraction for selected positions"));
 416     options->addOption(FileNameOption("olt").filetype(eftPlot).outputFile()
 417                            .store(&fnLifetime_).defaultBasename("lifetime")
 418                            .description("Lifetime histogram"));
 419
 420     selOpt_ = options->addOption(SelectionOption("select").storeVector(&sel_)
 421                                      .required().multiValue()
 422                                      .description("Selections to analyze"));
 423
 424     options->addOption(BooleanOption("norm").store(&bTotNorm_)
 425                            .description("Normalize by total number of positions with -os"));
 426     options->addOption(BooleanOption("cfnorm").store(&bFracNorm_)
 427                            .description("Normalize by covered fraction with -os"));
 428     const char *const cResNumberEnum[] = { "number", "index" };
 429     options->addOption(StringOption("resnr").store(&resNumberType_)
 430                            .enumValue(cResNumberEnum).defaultEnumIndex(0)
 431                            .description("Residue number output type with -oi and -on"));
 432     const char *const cPDBAtomsEnum[] = { "all", "maxsel", "selected" };
 433     options->addOption(StringOption("pdbatoms").store(&pdbAtoms_)
 434                            .enumValue(cPDBAtomsEnum).defaultEnumIndex(0)
 435                            .description("Atoms to write with -ofpdb"));
 436     options->addOption(BooleanOption("cumlt").store(&bCumulativeLifetimes_)
 437                            .description("Cumulate subintervals of longer intervals in -olt"));
 438 }
 439
 440 void
 441 Select::optionsFinished(TrajectoryAnalysisSettings *settings)
 442 {
 443     if (!fnPDB_.empty())
 444     {
 445         settings->setFlag(TrajectoryAnalysisSettings::efRequireTop);
 446         settings->setFlag(TrajectoryAnalysisSettings::efUseTopX);
 447     }
 448 }
 449
 450 void
 451 Select::initAnalysis(const TrajectoryAnalysisSettings &settings,
 452                      const TopologyInformation        &top)
 453 {
 454     bResInd_ = (resNumberType_ == "index");
 455
 456     for (SelectionList::iterator i = sel_.begin(); i != sel_.end(); ++i)
 457     {
 458         i->initCoveredFraction(CFRAC_SOLIDANGLE);
 459     }
 460
 461     // TODO: For large systems, a float may not have enough precision
 462     sdata_.setColumnCount(0, sel_.size());
 463     totsize_.reserve(sel_.size());
 464     for (size_t g = 0; g < sel_.size(); ++g)
 465     {
 466         totsize_.push_back(sel_[g].posCount());
 467     }
 468     if (!fnSize_.empty())
 469     {
 470         AnalysisDataPlotModulePointer plot(
 471                 new AnalysisDataPlotModule(settings.plotSettings()));
 472         plot->setFileName(fnSize_);
 473         plot->setTitle("Selection size");
 474         plot->setXAxisIsTime();
 475         plot->setYLabel("Number");
 476         for (size_t g = 0; g < sel_.size(); ++g)
 477         {
 478             plot->appendLegend(sel_[g].name());
 479         }
 480         sdata_.addModule(plot);
 481     }
 482
 483     cdata_.setColumnCount(0, sel_.size());
 484     if (!fnFrac_.empty())
 485     {
 486         AnalysisDataPlotModulePointer plot(
 487                 new AnalysisDataPlotModule(settings.plotSettings()));
 488         plot->setFileName(fnFrac_);
 489         plot->setTitle("Covered fraction");
 490         plot->setXAxisIsTime();
 491         plot->setYLabel("Fraction");
 492         plot->setYFormat(6, 4);
 493         for (size_t g = 0; g < sel_.size(); ++g)
 494         {
 495             plot->appendLegend(sel_[g].name());
 496         }
 497         cdata_.addModule(plot);
 498     }
 499
 500     // TODO: For large systems, a float may not have enough precision
 501     if (!fnIndex_.empty())
 502     {
 503         AnalysisDataPlotModulePointer plot(
 504                 new AnalysisDataPlotModule(settings.plotSettings()));
 505         plot->setFileName(fnIndex_);
 506         plot->setPlainOutput(true);
 507         plot->setYFormat(4, 0);
 508         idata_.addModule(plot);
 509     }
 510     if (!fnNdx_.empty())
 511     {
 512         boost::shared_ptr<IndexFileWriterModule> writer(new IndexFileWriterModule());
 513         writer->setFileName(fnNdx_);
 514         for (size_t g = 0; g < sel_.size(); ++g)
 515         {
 516             writer->addGroup(sel_[g].name(), sel_[g].isDynamic());
 517         }
 518         idata_.addModule(writer);
 519     }
 520
 521     mdata_.setDataSetCount(sel_.size());
 522     for (size_t g = 0; g < sel_.size(); ++g)
 523     {
 524         mdata_.setColumnCount(g, sel_[g].posCount());
 525     }
 526     lifetimeModule_->setCumulative(bCumulativeLifetimes_);
 527     if (!fnMask_.empty())
 528     {
 529         AnalysisDataPlotModulePointer plot(
 530                 new AnalysisDataPlotModule(settings.plotSettings()));
 531         plot->setFileName(fnMask_);
 532         plot->setTitle("Selection mask");
 533         plot->setXAxisIsTime();
 534         plot->setYLabel("Occupancy");
 535         plot->setYFormat(1, 0);
 536         // TODO: Add legend? (there can be massive amount of columns)
 537         mdata_.addModule(plot);
 538     }
 539     if (!fnOccupancy_.empty())
 540     {
 541         AnalysisDataPlotModulePointer plot(
 542                 new AnalysisDataPlotModule(settings.plotSettings()));
 543         plot->setFileName(fnOccupancy_);
 544         plot->setTitle("Fraction of time selection matches");
 545         plot->setXLabel("Selected position");
 546         plot->setYLabel("Occupied fraction");
 547         for (size_t g = 0; g < sel_.size(); ++g)
 548         {
 549             plot->appendLegend(sel_[g].name());
 550         }
 551         occupancyModule_->addModule(plot);
 552     }
 553     if (!fnLifetime_.empty())
 554     {
 555         AnalysisDataPlotModulePointer plot(
 556                 new AnalysisDataPlotModule(settings.plotSettings()));
 557         plot->setFileName(fnLifetime_);
 558         plot->setTitle("Lifetime histogram");
 559         plot->setXAxisIsTime();
 560         plot->setYLabel("Number of occurrences");
 561         for (size_t g = 0; g < sel_.size(); ++g)
 562         {
 563             plot->appendLegend(sel_[g].name());
 564         }
 565         lifetimeModule_->addModule(plot);
 566     }
 567
 568     top_ = &top;
 569 }
 570
 571
 572 void
 573 Select::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc * /* pbc */,
 574                      TrajectoryAnalysisModuleData *pdata)
 575 {
 576     AnalysisDataHandle   sdh = pdata->dataHandle(sdata_);
 577     AnalysisDataHandle   cdh = pdata->dataHandle(cdata_);
 578     AnalysisDataHandle   idh = pdata->dataHandle(idata_);
 579     AnalysisDataHandle   mdh = pdata->dataHandle(mdata_);
 580     const SelectionList &sel = pdata->parallelSelections(sel_);
 581     t_topology          *top = top_->topology();
 582
 583     sdh.startFrame(frnr, fr.time);
 584     for (size_t g = 0; g < sel.size(); ++g)
 585     {
 586         real normfac = bFracNorm_ ? 1.0 / sel[g].coveredFraction() : 1.0;
 587         if (bTotNorm_)
 588         {
 589             normfac /= totsize_[g];
 590         }
 591         sdh.setPoint(g, sel[g].posCount() * normfac);
 592     }
 593     sdh.finishFrame();
 594
 595     cdh.startFrame(frnr, fr.time);
 596     for (size_t g = 0; g < sel.size(); ++g)
 597     {
 598         cdh.setPoint(g, sel[g].coveredFraction());
 599     }
 600     cdh.finishFrame();
 601
 602     idh.startFrame(frnr, fr.time);
 603     for (size_t g = 0; g < sel.size(); ++g)
 604     {
 605         idh.setPoint(0, sel[g].posCount());
 606         idh.finishPointSet();
 607         for (int i = 0; i < sel[g].posCount(); ++i)
 608         {
 609             const SelectionPosition &p = sel[g].position(i);
 610             if (sel[g].type() == INDEX_RES && !bResInd_)
 611             {
 612                 idh.setPoint(1, top->atoms.resinfo[p.mappedId()].nr);
 613             }
 614             else
 615             {
 616                 idh.setPoint(1, p.mappedId() + 1);
 617             }
 618             idh.finishPointSet();
 619         }
 620     }
 621     idh.finishFrame();
 622
 623     mdh.startFrame(frnr, fr.time);
 624     for (size_t g = 0; g < sel.size(); ++g)
 625     {
 626         mdh.selectDataSet(g);
 627         for (int i = 0; i < totsize_[g]; ++i)
 628         {
 629             mdh.setPoint(i, 0);
 630         }
 631         for (int i = 0; i < sel[g].posCount(); ++i)
 632         {
 633             mdh.setPoint(sel[g].position(i).refId(), 1);
 634         }
 635     }
 636     mdh.finishFrame();
 637 }
 638
 639
 640 void
 641 Select::finishAnalysis(int /*nframes*/)
 642 {
 643 }
 644
 645
 646 void
 647 Select::writeOutput()
 648 {
 649     if (!fnPDB_.empty())
 650     {
 651         GMX_RELEASE_ASSERT(top_->hasTopology(),
 652                            "Topology should have been loaded or an error given earlier");
 653         t_atoms            atoms;
 654         atoms = top_->topology()->atoms;
 655         t_pdbinfo         *pdbinfo;
 656         snew(pdbinfo, atoms.nr);
 657         scoped_guard_sfree pdbinfoGuard(pdbinfo);
 658         if (atoms.pdbinfo != NULL)
 659         {
 660             std::memcpy(pdbinfo, atoms.pdbinfo, atoms.nr*sizeof(*pdbinfo));
 661         }
 662         atoms.pdbinfo = pdbinfo;
 663         for (int i = 0; i < atoms.nr; ++i)
 664         {
 665             pdbinfo[i].occup = 0.0;
 666         }
 667         for (size_t g = 0; g < sel_.size(); ++g)
 668         {
 669             for (int i = 0; i < sel_[g].posCount(); ++i)
 670             {
 671                 ConstArrayRef<int>                 atomIndices
 672                     = sel_[g].position(i).atomIndices();
 673                 ConstArrayRef<int>::const_iterator ai;
 674                 for (ai = atomIndices.begin(); ai != atomIndices.end(); ++ai)
 675                 {
 676                     pdbinfo[*ai].occup += occupancyModule_->average(g, i);
 677                 }
 678             }
 679         }
 680
 681         t_trxframe fr;
 682         clear_trxframe(&fr, TRUE);
 683         fr.bAtoms = TRUE;
 684         fr.atoms  = &atoms;
 685         fr.bX     = TRUE;
 686         fr.bBox   = TRUE;
 687         top_->getTopologyConf(&fr.x, fr.box);
 688
 689         if (pdbAtoms_ == "all")
 690         {
 691             t_trxstatus *status = open_trx(fnPDB_.c_str(), "w");
 692             write_trxframe(status, &fr, NULL);
 693             close_trx(status);
 694         }
 695         else if (pdbAtoms_ == "maxsel")
 696         {
 697             std::set<int> atomIndicesSet;
 698             for (size_t g = 0; g < sel_.size(); ++g)
 699             {
 700                 ConstArrayRef<int> atomIndices = sel_[g].atomIndices();
 701                 atomIndicesSet.insert(atomIndices.begin(), atomIndices.end());
 702             }
 703             std::vector<int>  allAtomIndices(atomIndicesSet.begin(),
 704                                              atomIndicesSet.end());
 705             t_trxstatus      *status = open_trx(fnPDB_.c_str(), "w");
 706             write_trxframe_indexed(status, &fr, allAtomIndices.size(),
 707                                    &allAtomIndices[0], NULL);
 708             close_trx(status);
 709         }
 710         else if (pdbAtoms_ == "selected")
 711         {
 712             std::vector<int> indices;
 713             for (int i = 0; i < atoms.nr; ++i)
 714             {
 715                 if (pdbinfo[i].occup > 0.0)
 716                 {
 717                     indices.push_back(i);
 718                 }
 719             }
 720             t_trxstatus *status = open_trx(fnPDB_.c_str(), "w");
 721             write_trxframe_indexed(status, &fr, indices.size(), &indices[0], NULL);
 722             close_trx(status);
 723         }
 724         else
 725         {
 726             GMX_RELEASE_ASSERT(false,
 727                                "Mismatch between -pdbatoms enum values and implementation");
 728         }
 729     }
 730 }
 731
 732 }       // namespace
 733
 734 const char SelectInfo::name[]             = "select";
 735 const char SelectInfo::shortDescription[] =
 736     "Print general information about selections";
 737
 738 TrajectoryAnalysisModulePointer SelectInfo::create()
 739 {
 740     return TrajectoryAnalysisModulePointer(new Select);
 741 }
 742
 743 } // namespace analysismodules
 744
 745 } // namespace gmx