src/gromacs/trajectoryanalysis/modules/trajectory.cpp

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