src/gromacs/trajectoryanalysis/modules/extract_cluster.cpp

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