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,2020, 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     output_env_init(&oenv_, getProgramContext(), settings_.timeUnit(), FALSE, XvgFormat::None, 0);
 212
 213     int frflags = settings_.frflags();
 214     frflags |= TRX_NEED_X;
 215
 216     snew(fr, 1);
 217
 218     if (hasTrajectory())
 219     {
 220         if (!read_first_frame(oenv_, &status_, trjfile_.c_str(), fr, frflags))
 221         {
 222             GMX_THROW(FileIOError("Could not read coordinates from trajectory"));
 223         }
 224         bTrajOpen_ = true;
 225
 226         if (topInfo_.hasTopology())
 227         {
 228             const int topologyAtomCount = topInfo_.mtop()->natoms;
 229             if (fr->natoms > topologyAtomCount)
 230             {
 231                 const std::string message =
 232                         formatString("Trajectory (%d atoms) does not match topology (%d atoms)",
 233                                      fr->natoms,
 234                                      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     setTrxFramePbcType(fr, topInfo_.pbcType());
 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(),
 285                 fr->natoms);
 286         GMX_THROW(InconsistentInputError(message));
 287     }
 288     fr->bIndex = TRUE;
 289     snew(fr->index, trajectoryGroup_.atomCount());
 290     std::copy(trajectoryGroup_.atomIndices().begin(), trajectoryGroup_.atomIndices().end(), fr->index);
 291 }
 292
 293 void TrajectoryAnalysisRunnerCommon::Impl::finishTrajectory()
 294 {
 295     if (bTrajOpen_)
 296     {
 297         close_trx(status_);
 298         bTrajOpen_ = false;
 299     }
 300     if (gpbc_ != nullptr)
 301     {
 302         gmx_rmpbc_done(gpbc_);
 303         gpbc_ = nullptr;
 304     }
 305 }
 306
 307 /*********************************************************************
 308  * TrajectoryAnalysisRunnerCommon
 309  */
 310
 311 TrajectoryAnalysisRunnerCommon::TrajectoryAnalysisRunnerCommon(TrajectoryAnalysisSettings* settings) :
 312     impl_(new Impl(settings))
 313 {
 314 }
 315
 316
 317 TrajectoryAnalysisRunnerCommon::~TrajectoryAnalysisRunnerCommon() {}
 318
 319
 320 ITopologyProvider* TrajectoryAnalysisRunnerCommon::topologyProvider()
 321 {
 322     return impl_.get();
 323 }
 324
 325
 326 void TrajectoryAnalysisRunnerCommon::initOptions(IOptionsContainer* options, TimeUnitBehavior* timeUnitBehavior)
 327 {
 328     TrajectoryAnalysisSettings& settings = impl_->settings_;
 329
 330     // Add common file name arguments.
 331     options->addOption(FileNameOption("f")
 332                                .filetype(eftTrajectory)
 333                                .inputFile()
 334                                .store(&impl_->trjfile_)
 335                                .defaultBasename("traj")
 336                                .description("Input trajectory or single configuration"));
 337     options->addOption(FileNameOption("s")
 338                                .filetype(eftTopology)
 339                                .inputFile()
 340                                .store(&impl_->topfile_)
 341                                .defaultBasename("topol")
 342                                .description("Input structure"));
 343
 344     // Add options for trajectory time control.
 345     options->addOption(DoubleOption("b")
 346                                .store(&impl_->startTime_)
 347                                .storeIsSet(&impl_->bStartTimeSet_)
 348                                .timeValue()
 349                                .description("First frame (%t) to read from trajectory"));
 350     options->addOption(DoubleOption("e")
 351                                .store(&impl_->endTime_)
 352                                .storeIsSet(&impl_->bEndTimeSet_)
 353                                .timeValue()
 354                                .description("Last frame (%t) to read from trajectory"));
 355     options->addOption(DoubleOption("dt")
 356                                .store(&impl_->deltaTime_)
 357                                .storeIsSet(&impl_->bDeltaTimeSet_)
 358                                .timeValue()
 359                                .description("Only use frame if t MOD dt == first time (%t)"));
 360
 361     // Add time unit option.
 362     timeUnitBehavior->setTimeUnitFromEnvironment();
 363     timeUnitBehavior->addTimeUnitOption(options, "tu");
 364     timeUnitBehavior->setTimeUnitStore(&impl_->settings_.impl_->timeUnit);
 365
 366     options->addOption(SelectionOption("fgroup")
 367                                .store(&impl_->trajectoryGroup_)
 368                                .onlySortedAtoms()
 369                                .onlyStatic()
 370                                .description("Atoms stored in the trajectory file "
 371                                             "(if not set, assume first N atoms)"));
 372
 373     // Add plot options.
 374     settings.impl_->plotSettings.initOptions(options);
 375
 376     // Add common options for trajectory processing.
 377     if (!settings.hasFlag(TrajectoryAnalysisSettings::efNoUserRmPBC))
 378     {
 379         options->addOption(
 380                 BooleanOption("rmpbc").store(&settings.impl_->bRmPBC).description("Make molecules whole for each frame"));
 381     }
 382     if (!settings.hasFlag(TrajectoryAnalysisSettings::efNoUserPBC))
 383     {
 384         options->addOption(
 385                 BooleanOption("pbc")
 386                         .store(&settings.impl_->bPBC)
 387                         .description("Use periodic boundary conditions for distance calculation"));
 388     }
 389 }
 390
 391
 392 void TrajectoryAnalysisRunnerCommon::optionsFinished()
 393 {
 394     if (impl_->trjfile_.empty() && impl_->topfile_.empty())
 395     {
 396         GMX_THROW(InconsistentInputError("No trajectory or topology provided, nothing to do!"));
 397     }
 398
 399     if (impl_->trajectoryGroup_.isValid() && impl_->trjfile_.empty())
 400     {
 401         GMX_THROW(
 402                 InconsistentInputError("-fgroup only makes sense together with a trajectory (-f)"));
 403     }
 404
 405     impl_->settings_.impl_->plotSettings.setTimeUnit(impl_->settings_.timeUnit());
 406
 407     if (impl_->bStartTimeSet_)
 408     {
 409         setTimeValue(TBEGIN, impl_->startTime_);
 410     }
 411     if (impl_->bEndTimeSet_)
 412     {
 413         setTimeValue(TEND, impl_->endTime_);
 414     }
 415     if (impl_->bDeltaTimeSet_)
 416     {
 417         setTimeValue(TDELTA, impl_->deltaTime_);
 418     }
 419 }
 420
 421
 422 void TrajectoryAnalysisRunnerCommon::initTopology()
 423 {
 424     const bool topologyRequired = impl_->settings_.hasFlag(TrajectoryAnalysisSettings::efRequireTop);
 425     impl_->initTopology(topologyRequired);
 426 }
 427
 428
 429 void TrajectoryAnalysisRunnerCommon::initFirstFrame()
 430 {
 431     impl_->initFirstFrame();
 432 }
 433
 434
 435 void TrajectoryAnalysisRunnerCommon::initFrameIndexGroup()
 436 {
 437     impl_->initFrameIndexGroup();
 438 }
 439
 440
 441 bool TrajectoryAnalysisRunnerCommon::readNextFrame()
 442 {
 443     bool bContinue = false;
 444     if (hasTrajectory())
 445     {
 446         bContinue = read_next_frame(impl_->oenv_, impl_->status_, impl_->fr);
 447     }
 448     if (!bContinue)
 449     {
 450         impl_->finishTrajectory();
 451     }
 452     return bContinue;
 453 }
 454
 455
 456 void TrajectoryAnalysisRunnerCommon::initFrame()
 457 {
 458     if (impl_->gpbc_ != nullptr)
 459     {
 460         gmx_rmpbc_trxfr(impl_->gpbc_, impl_->fr);
 461     }
 462 }
 463
 464
 465 bool TrajectoryAnalysisRunnerCommon::hasTrajectory() const
 466 {
 467     return impl_->hasTrajectory();
 468 }
 469
 470
 471 const TopologyInformation& TrajectoryAnalysisRunnerCommon::topologyInformation() const
 472 {
 473     return impl_->topInfo_;
 474 }
 475
 476
 477 t_trxframe& TrajectoryAnalysisRunnerCommon::frame() const
 478 {
 479     GMX_RELEASE_ASSERT(impl_->fr != nullptr, "Frame not available when accessed");
 480     return *impl_->fr;
 481 }
 482
 483 } // namespace gmx