src/gromacs/trajectoryanalysis/modules/trajectory.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2016,2018,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::Trajectory.
  38  *
  39  * \author Teemu Murtola <teemu.murtola@gmail.com>
  40  * \ingroup module_trajectoryanalysis
  41  */
  42 #include "gmxpre.h"
  43
  44 #include "trajectory.h"
  45
  46 #include <algorithm>
  47
  48 #include "gromacs/analysisdata/analysisdata.h"
  49 #include "gromacs/analysisdata/modules/plot.h"
  50 #include "gromacs/fileio/trxio.h"
  51 #include "gromacs/options/basicoptions.h"
  52 #include "gromacs/options/filenameoption.h"
  53 #include "gromacs/options/ioptionscontainer.h"
  54 #include "gromacs/selection/selection.h"
  55 #include "gromacs/selection/selectionoption.h"
  56 #include "gromacs/trajectory/trajectoryframe.h"
  57 #include "gromacs/trajectoryanalysis/analysissettings.h"
  58
  59 namespace gmx
  60 {
  61
  62 namespace analysismodules
  63 {
  64
  65 namespace
  66 {
  67
  68 /********************************************************************
  69  * Trajectory
  70  */
  71
  72 class Trajectory : public TrajectoryAnalysisModule
  73 {
  74 public:
  75     Trajectory();
  76
  77     void initOptions(IOptionsContainer* options, TrajectoryAnalysisSettings* settings) override;
  78     void optionsFinished(TrajectoryAnalysisSettings* settings) override;
  79     void initAnalysis(const TrajectoryAnalysisSettings& settings, const TopologyInformation& top) override;
  80
  81     void analyzeFrame(int frnr, const t_trxframe& fr, t_pbc* pbc, TrajectoryAnalysisModuleData* pdata) override;
  82
  83     void finishAnalysis(int nframes) override;
  84     void writeOutput() override;
  85
  86 private:
  87     SelectionList sel_;
  88
  89     std::string         fnX_;
  90     std::string         fnV_;
  91     std::string         fnF_;
  92     std::array<bool, 4> dimMask_;
  93     std::array<bool, 4> maskSet_;
  94
  95     AnalysisData xdata_;
  96     AnalysisData vdata_;
  97     AnalysisData fdata_;
  98 };
  99
 100 Trajectory::Trajectory() : dimMask_{ true, true, true, false }, maskSet_{}
 101 {
 102     registerAnalysisDataset(&xdata_, "x");
 103     registerAnalysisDataset(&vdata_, "v");
 104     registerAnalysisDataset(&fdata_, "f");
 105 }
 106
 107
 108 void Trajectory::initOptions(IOptionsContainer* options, TrajectoryAnalysisSettings* settings)
 109 {
 110     static const char* const desc[] = {
 111         "[THISMODULE] plots coordinates, velocities, and/or forces for",
 112         "provided selections. By default, the X, Y, and Z components for",
 113         "the requested vectors are plotted, but specifying one or more of",
 114         "[TT]-len[tt], [TT]-x[tt], [TT]-y[tt], and [TT]-z[tt] overrides this.",
 115         "",
 116         "For dynamic selections, currently the values are written out for",
 117         "all positions that the selection could select."
 118     };
 119
 120     settings->setHelpText(desc);
 121
 122     options->addOption(FileNameOption("ox")
 123                                .filetype(eftPlot)
 124                                .outputFile()
 125                                .store(&fnX_)
 126                                .defaultBasename("coord")
 127                                .description("Coordinates for each position as a function of time"));
 128     options->addOption(FileNameOption("ov")
 129                                .filetype(eftPlot)
 130                                .outputFile()
 131                                .store(&fnV_)
 132                                .defaultBasename("veloc")
 133                                .description("Velocities for each position as a function of time"));
 134     options->addOption(FileNameOption("of")
 135                                .filetype(eftPlot)
 136                                .outputFile()
 137                                .store(&fnF_)
 138                                .defaultBasename("force")
 139                                .description("Forces for each position as a function of time"));
 140
 141     options->addOption(
 142             SelectionOption("select").storeVector(&sel_).required().dynamicMask().multiValue().description(
 143                     "Selections to analyze"));
 144
 145     options->addOption(
 146             BooleanOption("x").store(&dimMask_[XX]).storeIsSet(&maskSet_[XX]).description("Plot X component"));
 147     options->addOption(
 148             BooleanOption("y").store(&dimMask_[YY]).storeIsSet(&maskSet_[YY]).description("Plot Y component"));
 149     options->addOption(
 150             BooleanOption("z").store(&dimMask_[ZZ]).storeIsSet(&maskSet_[ZZ]).description("Plot Z component"));
 151     options->addOption(
 152             BooleanOption("len").store(&dimMask_[DIM]).storeIsSet(&maskSet_[DIM]).description("Plot vector length"));
 153 }
 154
 155
 156 void Trajectory::optionsFinished(TrajectoryAnalysisSettings* settings)
 157 {
 158     int frameFlags = TRX_NEED_X;
 159     if (!fnV_.empty())
 160     {
 161         frameFlags |= TRX_READ_V;
 162     }
 163     if (!fnF_.empty())
 164     {
 165         frameFlags |= TRX_READ_F;
 166     }
 167     settings->setFrameFlags(frameFlags);
 168     if (std::count(std::begin(maskSet_), std::end(maskSet_), true) > 0)
 169     {
 170         for (int i = 0; i <= DIM; ++i)
 171         {
 172             if (!maskSet_[i])
 173             {
 174                 dimMask_[i] = false;
 175             }
 176         }
 177     }
 178 }
 179
 180
 181 void Trajectory::initAnalysis(const TrajectoryAnalysisSettings& settings, const TopologyInformation& /*top*/)
 182 {
 183     if (!fnX_.empty())
 184     {
 185         xdata_.setDataSetCount(sel_.size());
 186         for (size_t g = 0; g < sel_.size(); ++g)
 187         {
 188             xdata_.setColumnCount(g, 3 * sel_[g].posCount());
 189         }
 190         AnalysisDataVectorPlotModulePointer plot(new AnalysisDataVectorPlotModule(settings.plotSettings()));
 191         plot->setWriteMask(dimMask_.data());
 192         plot->setFileName(fnX_);
 193         plot->setTitle("Coordinates");
 194         plot->setXAxisIsTime();
 195         plot->setYLabel("Value [nm]");
 196         xdata_.addModule(plot);
 197     }
 198     if (!fnV_.empty())
 199     {
 200         vdata_.setDataSetCount(sel_.size());
 201         for (size_t g = 0; g < sel_.size(); ++g)
 202         {
 203             sel_[g].setEvaluateVelocities(true);
 204             vdata_.setColumnCount(g, 3 * sel_[g].posCount());
 205         }
 206         AnalysisDataVectorPlotModulePointer plot(new AnalysisDataVectorPlotModule(settings.plotSettings()));
 207         plot->setWriteMask(dimMask_.data());
 208         plot->setFileName(fnV_);
 209         plot->setTitle("Velocities");
 210         plot->setXAxisIsTime();
 211         plot->setYLabel("Value [nm/ps]");
 212         vdata_.addModule(plot);
 213     }
 214     if (!fnF_.empty())
 215     {
 216         fdata_.setDataSetCount(sel_.size());
 217         for (size_t g = 0; g < sel_.size(); ++g)
 218         {
 219             sel_[g].setEvaluateForces(true);
 220             fdata_.setColumnCount(g, 3 * sel_[g].posCount());
 221         }
 222         AnalysisDataVectorPlotModulePointer plot(new AnalysisDataVectorPlotModule(settings.plotSettings()));
 223         plot->setWriteMask(dimMask_.data());
 224         plot->setFileName(fnF_);
 225         plot->setTitle("Forces");
 226         plot->setXAxisIsTime();
 227         plot->setYLabel("Value [kJ mol\\S-1\\N nm\\S-1\\N]");
 228         fdata_.addModule(plot);
 229     }
 230 }
 231
 232
 233 //! Helper function for Trajectory::analyzeFrame
 234 template<typename T>
 235 void analyzeFrameImpl(int frnr, const t_trxframe& fr, AnalysisDataHandle* dh, const SelectionList& sel, T getField)
 236 {
 237     if (dh->isValid())
 238     {
 239         dh->startFrame(frnr, fr.time);
 240         for (size_t g = 0; g < sel.size(); ++g)
 241         {
 242             dh->selectDataSet(g);
 243             for (int i = 0; i < sel[g].posCount(); ++i)
 244             {
 245                 const SelectionPosition& pos = sel[g].position(i);
 246                 dh->setPoints(i * 3, 3, getField(pos), pos.selected());
 247             }
 248         }
 249         dh->finishFrame();
 250     }
 251 }
 252
 253 void Trajectory::analyzeFrame(int frnr, const t_trxframe& fr, t_pbc* /* pbc */, TrajectoryAnalysisModuleData* pdata)
 254 {
 255     AnalysisDataHandle   dh  = pdata->dataHandle(xdata_);
 256     const SelectionList& sel = TrajectoryAnalysisModuleData::parallelSelections(sel_);
 257     analyzeFrameImpl(frnr, fr, &dh, sel, [](const SelectionPosition& pos) { return pos.x(); });
 258     if (fr.bV)
 259     {
 260         analyzeFrameImpl(frnr, fr, &dh, sel, [](const SelectionPosition& pos) { return pos.v(); });
 261     }
 262     if (fr.bF)
 263     {
 264         analyzeFrameImpl(frnr, fr, &dh, sel, [](const SelectionPosition& pos) { return pos.f(); });
 265     }
 266 }
 267
 268
 269 void Trajectory::finishAnalysis(int /*nframes*/) {}
 270
 271
 272 void Trajectory::writeOutput() {}
 273
 274 } // namespace
 275
 276 const char TrajectoryInfo::name[] = "trajectory";
 277 const char TrajectoryInfo::shortDescription[] =
 278         "Print coordinates, velocities, and/or forces for selections";
 279
 280 TrajectoryAnalysisModulePointer TrajectoryInfo::create()
 281 {
 282     return TrajectoryAnalysisModulePointer(new Trajectory);
 283 }
 284
 285 } // namespace analysismodules
 286
 287 } // namespace gmx