src/gromacs/trajectoryanalysis/modules/extract_cluster.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2019, 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::ExtractCluster.
  38  *
  39  * \author Paul Bauer <paul.bauer.q@gmail.com>
  40  * \ingroup module_trajectoryanalysis
  41  */
  42
  43 #include "gmxpre.h"
  44
  45 #include "extract_cluster.h"
  46
  47 #include <algorithm>
  48
  49 #include "gromacs/coordinateio/coordinatefile.h"
  50 #include "gromacs/coordinateio/requirements.h"
  51 #include "gromacs/options/filenameoption.h"
  52 #include "gromacs/options/ioptionscontainer.h"
  53 #include "gromacs/selection/selectionoption.h"
  54 #include "gromacs/topology/index.h"
  55 #include "gromacs/trajectoryanalysis/analysissettings.h"
  56 #include "gromacs/trajectoryanalysis/topologyinformation.h"
  57 #include "gromacs/utility/path.h"
  58 #include "gromacs/utility/stringutil.h"
  59
  60 namespace gmx
  61 {
  62
  63 namespace analysismodules
  64 {
  65
  66 namespace
  67 {
  68
  69 /*
  70  * ExtractCluster
  71  */
  72
  73 class ExtractCluster : public TrajectoryAnalysisModule
  74 {
  75 public:
  76     ExtractCluster();
  77
  78     ~ExtractCluster() override;
  79
  80     void initOptions(IOptionsContainer* options, TrajectoryAnalysisSettings* settings) override;
  81     void optionsFinished(TrajectoryAnalysisSettings* settings) override;
  82     void initAnalysis(const TrajectoryAnalysisSettings& settings, const TopologyInformation& top) override;
  83     void analyzeFrame(int frameNumber, const t_trxframe& fr, t_pbc* pbc, TrajectoryAnalysisModuleData* pdata) override;
  84
  85     void finishAnalysis(int nframes) override;
  86     void writeOutput() override;
  87
  88 private:
  89     //! Storage of objects that handle output files.
  90     std::vector<TrajectoryFrameWriterPointer> writers_;
  91     //! Selection used for output.
  92     Selection sel_;
  93     //! Name for output file.
  94     std::string outputNamePrefix_;
  95     //! Name for index file.
  96     std::string indexFileName_;
  97     //! Storage of requirements for creating output files.
  98     OutputRequirementOptionDirector requirementsBuilder_;
  99     //! Stores the index information for the clusters. TODO refactor this!
 100     t_cluster_ndx* clusterIndex_ = nullptr;
 101 };
 102
 103 ExtractCluster::ExtractCluster() {}
 104
 105 ExtractCluster::~ExtractCluster()
 106 {
 107     if (clusterIndex_ != nullptr)
 108     {
 109         if (clusterIndex_->grpname != nullptr)
 110         {
 111             for (int i = 0; i < clusterIndex_->clust->nr; i++)
 112             {
 113                 sfree(clusterIndex_->grpname[i]);
 114             }
 115             sfree(clusterIndex_->grpname);
 116         }
 117         if (clusterIndex_->clust != nullptr)
 118         {
 119             sfree(clusterIndex_->inv_clust);
 120             done_blocka(clusterIndex_->clust);
 121             sfree(clusterIndex_->clust);
 122         }
 123         sfree(clusterIndex_);
 124     }
 125 }
 126
 127
 128 void ExtractCluster::initOptions(IOptionsContainer* options, TrajectoryAnalysisSettings* settings)
 129 {
 130     static const char* const desc[] = {
 131         "[THISMODULE] can be used to extract trajectory frames that correspond to clusters ",
 132         "obtained from running gmx cluster with the -clndx option.",
 133         "The module supports writing all GROMACS supported trajectory file formats.",
 134         "",
 135         "Included is also a selection of possible options to change additional information.",
 136         "",
 137         "It is possible to write only a selection of atoms to the output trajectory",
 138         "files for each cluster.",
 139     };
 140
 141     options->addOption(FileNameOption("clusters")
 142                                .filetype(eftIndex)
 143                                .inputFile()
 144                                .required()
 145                                .store(&indexFileName_)
 146                                .defaultBasename("cluster")
 147                                .description("Name of index file containing frame indices for each "
 148                                             "cluster, obtained from gmx cluster -clndx."));
 149
 150     options->addOption(SelectionOption("select").store(&sel_).onlyAtoms().description(
 151             "Selection of atoms to write to the file"));
 152
 153     options->addOption(FileNameOption("o")
 154                                .filetype(eftTrajectory)
 155                                .outputFile()
 156                                .store(&outputNamePrefix_)
 157                                .defaultBasename("trajout")
 158                                .required()
 159                                .description("Prefix for the name of the trajectory file written "
 160                                             "for each cluster."));
 161
 162     requirementsBuilder_.initOptions(options);
 163
 164     settings->setHelpText(desc);
 165 }
 166
 167 void ExtractCluster::optionsFinished(TrajectoryAnalysisSettings* settings)
 168 {
 169     int frameFlags = TRX_NEED_X;
 170
 171     frameFlags |= TRX_READ_V | TRX_READ_F;
 172
 173     settings->setFrameFlags(frameFlags);
 174
 175     clusterIndex_ = cluster_index(nullptr, indexFileName_.c_str());
 176 }
 177
 178
 179 void ExtractCluster::initAnalysis(const TrajectoryAnalysisSettings& /*settings*/,
 180                                   const TopologyInformation& top)
 181 {
 182     int numberOfClusters = clusterIndex_->clust->nr;
 183     for (int i = 0; i < numberOfClusters; i++)
 184     {
 185         std::string outputName = Path::concatenateBeforeExtension(
 186                 outputNamePrefix_, formatString("_%s", clusterIndex_->grpname[i]));
 187         writers_.emplace_back(createTrajectoryFrameWriter(top.mtop(), sel_, outputName,
 188                                                           top.hasTopology() ? top.copyAtoms() : nullptr,
 189                                                           requirementsBuilder_.process()));
 190     }
 191 }
 192
 193 void ExtractCluster::analyzeFrame(int               frameNumber,
 194                                   const t_trxframe& frame,
 195                                   t_pbc* /* pbc */,
 196                                   TrajectoryAnalysisModuleData* /*pdata*/)
 197 {
 198     // modify frame to write out correct number of coords
 199     // and actually write out
 200     int clusterToWriteTo = clusterIndex_->inv_clust[frameNumber];
 201     // Check for valid entry in cluster list, otherwise skip frame.
 202     if (clusterToWriteTo != -1 && clusterToWriteTo < clusterIndex_->clust->nr)
 203     {
 204         writers_[clusterToWriteTo]->prepareAndWriteFrame(frameNumber, frame);
 205     }
 206     else
 207     {
 208         printf("Frame %d was not found in any cluster!", frameNumber);
 209     }
 210 }
 211
 212 void ExtractCluster::finishAnalysis(int /*nframes*/) {}
 213
 214
 215 void ExtractCluster::writeOutput() {}
 216
 217 } // namespace
 218
 219 const char ExtractClusterInfo::name[] = "extract-cluster";
 220 const char ExtractClusterInfo::shortDescription[] =
 221         "Allows extracting frames corresponding to clusters from trajectory";
 222
 223 TrajectoryAnalysisModulePointer ExtractClusterInfo::create()
 224 {
 225     return TrajectoryAnalysisModulePointer(new ExtractCluster);
 226 }
 227
 228 } // namespace analysismodules
 229
 230 } // namespace gmx