src/gromacs/trajectoryanalysis/runnercommon.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2010-2018, The GROMACS development team.
   5  * Copyright (c) 2019, by the GROMACS development team, led by
   6  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   7  * and including many others, as listed in the AUTHORS file in the
   8  * top-level source directory and at http://www.gromacs.org.
   9  *
  10  * GROMACS is free software; you can redistribute it and/or
  11  * modify it under the terms of the GNU Lesser General Public License
  12  * as published by the Free Software Foundation; either version 2.1
  13  * of the License, or (at your option) any later version.
  14  *
  15  * GROMACS is distributed in the hope that it will be useful,
  16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  18  * Lesser General Public License for more details.
  19  *
  20  * You should have received a copy of the GNU Lesser General Public
  21  * License along with GROMACS; if not, see
  22  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  23  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  24  *
  25  * If you want to redistribute modifications to GROMACS, please
  26  * consider that scientific software is very special. Version
  27  * control is crucial - bugs must be traceable. We will be happy to
  28  * consider code for inclusion in the official distribution, but
  29  * derived work must not be called official GROMACS. Details are found
  30  * in the README & COPYING files - if they are missing, get the
  31  * official version at http://www.gromacs.org.
  32  *
  33  * To help us fund GROMACS development, we humbly ask that you cite
  34  * the research papers on the package. Check out http://www.gromacs.org.
  35  */
  36 /*! \internal \file
  37  * \brief
  38  * Implements gmx::TrajectoryAnalysisRunnerCommon.
  39  *
  40  * \author Teemu Murtola <teemu.murtola@gmail.com>
  41  * \ingroup module_trajectoryanalysis
  42  */
  43 #include "gmxpre.h"
  44
  45 #include "runnercommon.h"
  46
  47 #include <cstring>
  48
  49 #include <algorithm>
  50 #include <string>
  51
  52 #include "gromacs/fileio/oenv.h"
  53 #include "gromacs/fileio/timecontrol.h"
  54 #include "gromacs/fileio/trxio.h"
  55 #include "gromacs/math/vec.h"
  56 #include "gromacs/options/basicoptions.h"
  57 #include "gromacs/options/filenameoption.h"
  58 #include "gromacs/options/ioptionscontainer.h"
  59 #include "gromacs/pbcutil/rmpbc.h"
  60 #include "gromacs/selection/selection.h"
  61 #include "gromacs/selection/selectioncollection.h"
  62 #include "gromacs/selection/selectionoption.h"
  63 #include "gromacs/selection/selectionoptionbehavior.h"
  64 #include "gromacs/topology/topology.h"
  65 #include "gromacs/trajectory/trajectoryframe.h"
  66 #include "gromacs/trajectoryanalysis/analysissettings.h"
  67 #include "gromacs/trajectoryanalysis/topologyinformation.h"
  68 #include "gromacs/utility/cstringutil.h"
  69 #include "gromacs/utility/exceptions.h"
  70 #include "gromacs/utility/gmxassert.h"
  71 #include "gromacs/utility/programcontext.h"
  72 #include "gromacs/utility/smalloc.h"
  73 #include "gromacs/utility/stringutil.h"
  74
  75 #include "analysissettings_impl.h"
  76
  77 namespace gmx
  78 {
  79
  80 class TrajectoryAnalysisRunnerCommon::Impl : public ITopologyProvider
  81 {
  82 public:
  83     explicit Impl(TrajectoryAnalysisSettings* settings);
  84     ~Impl() override;
  85
  86     bool hasTrajectory() const { return !trjfile_.empty(); }
  87
  88     void initTopology(bool required);
  89     void initFirstFrame();
  90     void initFrameIndexGroup();
  91     void finishTrajectory();
  92
  93     // From ITopologyProvider
  94     gmx_mtop_t* getTopology(bool required) override
  95     {
  96         initTopology(required);
  97         return topInfo_.mtop_.get();
  98     }
  99     int getAtomCount() override
 100     {
 101         if (!topInfo_.hasTopology())
 102         {
 103             if (trajectoryGroup_.isValid())
 104             {
 105                 GMX_THROW(InconsistentInputError(
 106                         "-fgroup is only supported when -s is also specified"));
 107             }
 108             // Read the first frame if we don't know the maximum number of
 109             // atoms otherwise.
 110             initFirstFrame();
 111             return fr->natoms;
 112         }
 113         return -1;
 114     }
 115
 116     TrajectoryAnalysisSettings& settings_;
 117     TopologyInformation         topInfo_;
 118
 119     //! Name of the trajectory file (empty if not provided).
 120     std::string trjfile_;
 121     //! Name of the topology file (empty if no topology provided).
 122     std::string topfile_;
 123     Selection   trajectoryGroup_;
 124     double      startTime_;
 125     double      endTime_;
 126     double      deltaTime_;
 127     bool        bStartTimeSet_;
 128     bool        bEndTimeSet_;
 129     bool        bDeltaTimeSet_;
 130
 131     bool bTrajOpen_;
 132     //! The current frame, or \p NULL if no frame loaded yet.
 133     t_trxframe* fr;
 134     gmx_rmpbc_t gpbc_;
 135     //! Used to store the status variable from read_first_frame().
 136     t_trxstatus*      status_;
 137     gmx_output_env_t* oenv_;
 138 };
 139
 140
 141 TrajectoryAnalysisRunnerCommon::Impl::Impl(TrajectoryAnalysisSettings* settings) :
 142     settings_(*settings),
 143     startTime_(0.0),
 144     endTime_(0.0),
 145     deltaTime_(0.0),
 146     bStartTimeSet_(false),
 147     bEndTimeSet_(false),
 148     bDeltaTimeSet_(false),
 149     bTrajOpen_(false),
 150     fr(nullptr),
 151     gpbc_(nullptr),
 152     status_(nullptr),
 153     oenv_(nullptr)
 154 {
 155 }
 156
 157
 158 TrajectoryAnalysisRunnerCommon::Impl::~Impl()
 159 {
 160     finishTrajectory();
 161     if (fr != nullptr)
 162     {
 163         // There doesn't seem to be a function for freeing frame data
 164         sfree(fr->x);
 165         sfree(fr->v);
 166         sfree(fr->f);
 167         sfree(fr->index);
 168         sfree(fr);
 169     }
 170     if (oenv_ != nullptr)
 171     {
 172         output_env_done(oenv_);
 173     }
 174 }
 175
 176 void TrajectoryAnalysisRunnerCommon::Impl::initTopology(bool required)
 177 {
 178     // Return immediately if the topology has already been loaded.
 179     if (topInfo_.hasTopology())
 180     {
 181         return;
 182     }
 183
 184     if (required && topfile_.empty())
 185     {
 186         GMX_THROW(InconsistentInputError("No topology provided, but one is required for analysis"));
 187     }
 188
 189     // Load the topology if requested.
 190     if (!topfile_.empty())
 191     {
 192         topInfo_.fillFromInputFile(topfile_);
 193         if (hasTrajectory() && !settings_.hasFlag(TrajectoryAnalysisSettings::efUseTopX))
 194         {
 195             topInfo_.xtop_.clear();
 196         }
 197         if (hasTrajectory() && !settings_.hasFlag(TrajectoryAnalysisSettings::efUseTopV))
 198         {
 199             topInfo_.vtop_.clear();
 200         }
 201     }
 202 }
 203
 204 void TrajectoryAnalysisRunnerCommon::Impl::initFirstFrame()
 205 {
 206     // Return if we have already initialized the trajectory.
 207     if (fr != nullptr)
 208     {
 209         return;
 210     }
 211     time_unit_t time_unit = static_cast<time_unit_t>(settings_.timeUnit() + 1); // NOLINT(bugprone-misplaced-widening-cast)
 212     output_env_init(&oenv_, getProgramContext(), time_unit, FALSE, exvgNONE, 0);
 213
 214     int frflags = settings_.frflags();
 215     frflags |= TRX_NEED_X;
 216
 217     snew(fr, 1);
 218
 219     if (hasTrajectory())
 220     {
 221         if (!read_first_frame(oenv_, &status_, trjfile_.c_str(), fr, frflags))
 222         {
 223             GMX_THROW(FileIOError("Could not read coordinates from trajectory"));
 224         }
 225         bTrajOpen_ = true;
 226
 227         if (topInfo_.hasTopology())
 228         {
 229             const int topologyAtomCount = topInfo_.mtop()->natoms;
 230             if (fr->natoms > topologyAtomCount)
 231             {
 232                 const std::string message =
 233                         formatString("Trajectory (%d atoms) does not match topology (%d atoms)",
 234                                      fr->natoms, topologyAtomCount);
 235                 GMX_THROW(InconsistentInputError(message));
 236             }
 237         }
 238     }
 239     else
 240     {
 241         // Prepare a frame from topology information.
 242         if (frflags & (TRX_NEED_F))
 243         {
 244             GMX_THROW(InvalidInputError("Forces cannot be read from a topology"));
 245         }
 246         fr->natoms = topInfo_.mtop()->natoms;
 247         fr->bX     = TRUE;
 248         snew(fr->x, fr->natoms);
 249         memcpy(fr->x, topInfo_.xtop_.data(), sizeof(*fr->x) * fr->natoms);
 250         if (frflags & (TRX_NEED_V))
 251         {
 252             if (topInfo_.vtop_.empty())
 253             {
 254                 GMX_THROW(InvalidInputError(
 255                         "Velocities were required, but could not be read from the topology file"));
 256             }
 257             fr->bV = TRUE;
 258             snew(fr->v, fr->natoms);
 259             memcpy(fr->v, topInfo_.vtop_.data(), sizeof(*fr->v) * fr->natoms);
 260         }
 261         fr->bBox = TRUE;
 262         copy_mat(topInfo_.boxtop_, fr->box);
 263     }
 264
 265     set_trxframe_ePBC(fr, topInfo_.ePBC());
 266     if (topInfo_.hasTopology() && settings_.hasRmPBC())
 267     {
 268         gpbc_ = gmx_rmpbc_init(topInfo_);
 269     }
 270 }
 271
 272 void TrajectoryAnalysisRunnerCommon::Impl::initFrameIndexGroup()
 273 {
 274     if (!trajectoryGroup_.isValid())
 275     {
 276         return;
 277     }
 278     GMX_RELEASE_ASSERT(bTrajOpen_, "Trajectory index only makes sense with a real trajectory");
 279     if (trajectoryGroup_.atomCount() != fr->natoms)
 280     {
 281         const std::string message = formatString(
 282                 "Selection specified with -fgroup has %d atoms, but "
 283                 "the trajectory (-f) has %d atoms.",
 284                 trajectoryGroup_.atomCount(), fr->natoms);
 285         GMX_THROW(InconsistentInputError(message));
 286     }
 287     fr->bIndex = TRUE;
 288     snew(fr->index, trajectoryGroup_.atomCount());
 289     std::copy(trajectoryGroup_.atomIndices().begin(), trajectoryGroup_.atomIndices().end(), fr->index);
 290 }
 291
 292 void TrajectoryAnalysisRunnerCommon::Impl::finishTrajectory()
 293 {
 294     if (bTrajOpen_)
 295     {
 296         close_trx(status_);
 297         bTrajOpen_ = false;
 298     }
 299     if (gpbc_ != nullptr)
 300     {
 301         gmx_rmpbc_done(gpbc_);
 302         gpbc_ = nullptr;
 303     }
 304 }
 305
 306 /*********************************************************************
 307  * TrajectoryAnalysisRunnerCommon
 308  */
 309
 310 TrajectoryAnalysisRunnerCommon::TrajectoryAnalysisRunnerCommon(TrajectoryAnalysisSettings* settings) :
 311     impl_(new Impl(settings))
 312 {
 313 }
 314
 315
 316 TrajectoryAnalysisRunnerCommon::~TrajectoryAnalysisRunnerCommon() {}
 317
 318
 319 ITopologyProvider* TrajectoryAnalysisRunnerCommon::topologyProvider()
 320 {
 321     return impl_.get();
 322 }
 323
 324
 325 void TrajectoryAnalysisRunnerCommon::initOptions(IOptionsContainer* options, TimeUnitBehavior* timeUnitBehavior)
 326 {
 327     TrajectoryAnalysisSettings& settings = impl_->settings_;
 328
 329     // Add common file name arguments.
 330     options->addOption(FileNameOption("f")
 331                                .filetype(eftTrajectory)
 332                                .inputFile()
 333                                .store(&impl_->trjfile_)
 334                                .defaultBasename("traj")
 335                                .description("Input trajectory or single configuration"));
 336     options->addOption(FileNameOption("s")
 337                                .filetype(eftTopology)
 338                                .inputFile()
 339                                .store(&impl_->topfile_)
 340                                .defaultBasename("topol")
 341                                .description("Input structure"));
 342
 343     // Add options for trajectory time control.
 344     options->addOption(DoubleOption("b")
 345                                .store(&impl_->startTime_)
 346                                .storeIsSet(&impl_->bStartTimeSet_)
 347                                .timeValue()
 348                                .description("First frame (%t) to read from trajectory"));
 349     options->addOption(DoubleOption("e")
 350                                .store(&impl_->endTime_)
 351                                .storeIsSet(&impl_->bEndTimeSet_)
 352                                .timeValue()
 353                                .description("Last frame (%t) to read from trajectory"));
 354     options->addOption(DoubleOption("dt")
 355                                .store(&impl_->deltaTime_)
 356                                .storeIsSet(&impl_->bDeltaTimeSet_)
 357                                .timeValue()
 358                                .description("Only use frame if t MOD dt == first time (%t)"));
 359
 360     // Add time unit option.
 361     timeUnitBehavior->setTimeUnitFromEnvironment();
 362     timeUnitBehavior->addTimeUnitOption(options, "tu");
 363     timeUnitBehavior->setTimeUnitStore(&impl_->settings_.impl_->timeUnit);
 364
 365     options->addOption(SelectionOption("fgroup")
 366                                .store(&impl_->trajectoryGroup_)
 367                                .onlySortedAtoms()
 368                                .onlyStatic()
 369                                .description("Atoms stored in the trajectory file "
 370                                             "(if not set, assume first N atoms)"));
 371
 372     // Add plot options.
 373     settings.impl_->plotSettings.initOptions(options);
 374
 375     // Add common options for trajectory processing.
 376     if (!settings.hasFlag(TrajectoryAnalysisSettings::efNoUserRmPBC))
 377     {
 378         options->addOption(
 379                 BooleanOption("rmpbc").store(&settings.impl_->bRmPBC).description("Make molecules whole for each frame"));
 380     }
 381     if (!settings.hasFlag(TrajectoryAnalysisSettings::efNoUserPBC))
 382     {
 383         options->addOption(
 384                 BooleanOption("pbc")
 385                         .store(&settings.impl_->bPBC)
 386                         .description("Use periodic boundary conditions for distance calculation"));
 387     }
 388 }
 389
 390
 391 void TrajectoryAnalysisRunnerCommon::optionsFinished()
 392 {
 393     if (impl_->trjfile_.empty() && impl_->topfile_.empty())
 394     {
 395         GMX_THROW(InconsistentInputError("No trajectory or topology provided, nothing to do!"));
 396     }
 397
 398     if (impl_->trajectoryGroup_.isValid() && impl_->trjfile_.empty())
 399     {
 400         GMX_THROW(
 401                 InconsistentInputError("-fgroup only makes sense together with a trajectory (-f)"));
 402     }
 403
 404     impl_->settings_.impl_->plotSettings.setTimeUnit(impl_->settings_.timeUnit());
 405
 406     if (impl_->bStartTimeSet_)
 407     {
 408         setTimeValue(TBEGIN, impl_->startTime_);
 409     }
 410     if (impl_->bEndTimeSet_)
 411     {
 412         setTimeValue(TEND, impl_->endTime_);
 413     }
 414     if (impl_->bDeltaTimeSet_)
 415     {
 416         setTimeValue(TDELTA, impl_->deltaTime_);
 417     }
 418 }
 419
 420
 421 void TrajectoryAnalysisRunnerCommon::initTopology()
 422 {
 423     const bool topologyRequired = impl_->settings_.hasFlag(TrajectoryAnalysisSettings::efRequireTop);
 424     impl_->initTopology(topologyRequired);
 425 }
 426
 427
 428 void TrajectoryAnalysisRunnerCommon::initFirstFrame()
 429 {
 430     impl_->initFirstFrame();
 431 }
 432
 433
 434 void TrajectoryAnalysisRunnerCommon::initFrameIndexGroup()
 435 {
 436     impl_->initFrameIndexGroup();
 437 }
 438
 439
 440 bool TrajectoryAnalysisRunnerCommon::readNextFrame()
 441 {
 442     bool bContinue = false;
 443     if (hasTrajectory())
 444     {
 445         bContinue = read_next_frame(impl_->oenv_, impl_->status_, impl_->fr);
 446     }
 447     if (!bContinue)
 448     {
 449         impl_->finishTrajectory();
 450     }
 451     return bContinue;
 452 }
 453
 454
 455 void TrajectoryAnalysisRunnerCommon::initFrame()
 456 {
 457     if (impl_->gpbc_ != nullptr)
 458     {
 459         gmx_rmpbc_trxfr(impl_->gpbc_, impl_->fr);
 460     }
 461 }
 462
 463
 464 bool TrajectoryAnalysisRunnerCommon::hasTrajectory() const
 465 {
 466     return impl_->hasTrajectory();
 467 }
 468
 469
 470 const TopologyInformation& TrajectoryAnalysisRunnerCommon::topologyInformation() const
 471 {
 472     return impl_->topInfo_;
 473 }
 474
 475
 476 t_trxframe& TrajectoryAnalysisRunnerCommon::frame() const
 477 {
 478     GMX_RELEASE_ASSERT(impl_->fr != nullptr, "Frame not available when accessed");
 479     return *impl_->fr;
 480 }
 481
 482 } // namespace gmx