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, topologyAtomCount);
 234                 GMX_THROW(InconsistentInputError(message));
 235             }
 236         }
 237     }
 238     else
 239     {
 240         // Prepare a frame from topology information.
 241         if (frflags & (TRX_NEED_F))
 242         {
 243             GMX_THROW(InvalidInputError("Forces cannot be read from a topology"));
 244         }
 245         fr->natoms = topInfo_.mtop()->natoms;
 246         fr->bX     = TRUE;
 247         snew(fr->x, fr->natoms);
 248         memcpy(fr->x, topInfo_.xtop_.data(), sizeof(*fr->x) * fr->natoms);
 249         if (frflags & (TRX_NEED_V))
 250         {
 251             if (topInfo_.vtop_.empty())
 252             {
 253                 GMX_THROW(InvalidInputError(
 254                         "Velocities were required, but could not be read from the topology file"));
 255             }
 256             fr->bV = TRUE;
 257             snew(fr->v, fr->natoms);
 258             memcpy(fr->v, topInfo_.vtop_.data(), sizeof(*fr->v) * fr->natoms);
 259         }
 260         fr->bBox = TRUE;
 261         copy_mat(topInfo_.boxtop_, fr->box);
 262     }
 263
 264     setTrxFramePbcType(fr, topInfo_.pbcType());
 265     if (topInfo_.hasTopology() && settings_.hasRmPBC())
 266     {
 267         gpbc_ = gmx_rmpbc_init(topInfo_);
 268     }
 269 }
 270
 271 void TrajectoryAnalysisRunnerCommon::Impl::initFrameIndexGroup()
 272 {
 273     if (!trajectoryGroup_.isValid())
 274     {
 275         return;
 276     }
 277     GMX_RELEASE_ASSERT(bTrajOpen_, "Trajectory index only makes sense with a real trajectory");
 278     if (trajectoryGroup_.atomCount() != fr->natoms)
 279     {
 280         const std::string message = formatString(
 281                 "Selection specified with -fgroup has %d atoms, but "
 282                 "the trajectory (-f) has %d atoms.",
 283                 trajectoryGroup_.atomCount(), fr->natoms);
 284         GMX_THROW(InconsistentInputError(message));
 285     }
 286     fr->bIndex = TRUE;
 287     snew(fr->index, trajectoryGroup_.atomCount());
 288     std::copy(trajectoryGroup_.atomIndices().begin(), trajectoryGroup_.atomIndices().end(), fr->index);
 289 }
 290
 291 void TrajectoryAnalysisRunnerCommon::Impl::finishTrajectory()
 292 {
 293     if (bTrajOpen_)
 294     {
 295         close_trx(status_);
 296         bTrajOpen_ = false;
 297     }
 298     if (gpbc_ != nullptr)
 299     {
 300         gmx_rmpbc_done(gpbc_);
 301         gpbc_ = nullptr;
 302     }
 303 }
 304
 305 /*********************************************************************
 306  * TrajectoryAnalysisRunnerCommon
 307  */
 308
 309 TrajectoryAnalysisRunnerCommon::TrajectoryAnalysisRunnerCommon(TrajectoryAnalysisSettings* settings) :
 310     impl_(new Impl(settings))
 311 {
 312 }
 313
 314
 315 TrajectoryAnalysisRunnerCommon::~TrajectoryAnalysisRunnerCommon() {}
 316
 317
 318 ITopologyProvider* TrajectoryAnalysisRunnerCommon::topologyProvider()
 319 {
 320     return impl_.get();
 321 }
 322
 323
 324 void TrajectoryAnalysisRunnerCommon::initOptions(IOptionsContainer* options, TimeUnitBehavior* timeUnitBehavior)
 325 {
 326     TrajectoryAnalysisSettings& settings = impl_->settings_;
 327
 328     // Add common file name arguments.
 329     options->addOption(FileNameOption("f")
 330                                .filetype(eftTrajectory)
 331                                .inputFile()
 332                                .store(&impl_->trjfile_)
 333                                .defaultBasename("traj")
 334                                .description("Input trajectory or single configuration"));
 335     options->addOption(FileNameOption("s")
 336                                .filetype(eftTopology)
 337                                .inputFile()
 338                                .store(&impl_->topfile_)
 339                                .defaultBasename("topol")
 340                                .description("Input structure"));
 341
 342     // Add options for trajectory time control.
 343     options->addOption(DoubleOption("b")
 344                                .store(&impl_->startTime_)
 345                                .storeIsSet(&impl_->bStartTimeSet_)
 346                                .timeValue()
 347                                .description("First frame (%t) to read from trajectory"));
 348     options->addOption(DoubleOption("e")
 349                                .store(&impl_->endTime_)
 350                                .storeIsSet(&impl_->bEndTimeSet_)
 351                                .timeValue()
 352                                .description("Last frame (%t) to read from trajectory"));
 353     options->addOption(DoubleOption("dt")
 354                                .store(&impl_->deltaTime_)
 355                                .storeIsSet(&impl_->bDeltaTimeSet_)
 356                                .timeValue()
 357                                .description("Only use frame if t MOD dt == first time (%t)"));
 358
 359     // Add time unit option.
 360     timeUnitBehavior->setTimeUnitFromEnvironment();
 361     timeUnitBehavior->addTimeUnitOption(options, "tu");
 362     timeUnitBehavior->setTimeUnitStore(&impl_->settings_.impl_->timeUnit);
 363
 364     options->addOption(SelectionOption("fgroup")
 365                                .store(&impl_->trajectoryGroup_)
 366                                .onlySortedAtoms()
 367                                .onlyStatic()
 368                                .description("Atoms stored in the trajectory file "
 369                                             "(if not set, assume first N atoms)"));
 370
 371     // Add plot options.
 372     settings.impl_->plotSettings.initOptions(options);
 373
 374     // Add common options for trajectory processing.
 375     if (!settings.hasFlag(TrajectoryAnalysisSettings::efNoUserRmPBC))
 376     {
 377         options->addOption(
 378                 BooleanOption("rmpbc").store(&settings.impl_->bRmPBC).description("Make molecules whole for each frame"));
 379     }
 380     if (!settings.hasFlag(TrajectoryAnalysisSettings::efNoUserPBC))
 381     {
 382         options->addOption(
 383                 BooleanOption("pbc")
 384                         .store(&settings.impl_->bPBC)
 385                         .description("Use periodic boundary conditions for distance calculation"));
 386     }
 387 }
 388
 389
 390 void TrajectoryAnalysisRunnerCommon::optionsFinished()
 391 {
 392     if (impl_->trjfile_.empty() && impl_->topfile_.empty())
 393     {
 394         GMX_THROW(InconsistentInputError("No trajectory or topology provided, nothing to do!"));
 395     }
 396
 397     if (impl_->trajectoryGroup_.isValid() && impl_->trjfile_.empty())
 398     {
 399         GMX_THROW(
 400                 InconsistentInputError("-fgroup only makes sense together with a trajectory (-f)"));
 401     }
 402
 403     impl_->settings_.impl_->plotSettings.setTimeUnit(impl_->settings_.timeUnit());
 404
 405     if (impl_->bStartTimeSet_)
 406     {
 407         setTimeValue(TBEGIN, impl_->startTime_);
 408     }
 409     if (impl_->bEndTimeSet_)
 410     {
 411         setTimeValue(TEND, impl_->endTime_);
 412     }
 413     if (impl_->bDeltaTimeSet_)
 414     {
 415         setTimeValue(TDELTA, impl_->deltaTime_);
 416     }
 417 }
 418
 419
 420 void TrajectoryAnalysisRunnerCommon::initTopology()
 421 {
 422     const bool topologyRequired = impl_->settings_.hasFlag(TrajectoryAnalysisSettings::efRequireTop);
 423     impl_->initTopology(topologyRequired);
 424 }
 425
 426
 427 void TrajectoryAnalysisRunnerCommon::initFirstFrame()
 428 {
 429     impl_->initFirstFrame();
 430 }
 431
 432
 433 void TrajectoryAnalysisRunnerCommon::initFrameIndexGroup()
 434 {
 435     impl_->initFrameIndexGroup();
 436 }
 437
 438
 439 bool TrajectoryAnalysisRunnerCommon::readNextFrame()
 440 {
 441     bool bContinue = false;
 442     if (hasTrajectory())
 443     {
 444         bContinue = read_next_frame(impl_->oenv_, impl_->status_, impl_->fr);
 445     }
 446     if (!bContinue)
 447     {
 448         impl_->finishTrajectory();
 449     }
 450     return bContinue;
 451 }
 452
 453
 454 void TrajectoryAnalysisRunnerCommon::initFrame()
 455 {
 456     if (impl_->gpbc_ != nullptr)
 457     {
 458         gmx_rmpbc_trxfr(impl_->gpbc_, impl_->fr);
 459     }
 460 }
 461
 462
 463 bool TrajectoryAnalysisRunnerCommon::hasTrajectory() const
 464 {
 465     return impl_->hasTrajectory();
 466 }
 467
 468
 469 const TopologyInformation& TrajectoryAnalysisRunnerCommon::topologyInformation() const
 470 {
 471     return impl_->topInfo_;
 472 }
 473
 474
 475 t_trxframe& TrajectoryAnalysisRunnerCommon::frame() const
 476 {
 477     GMX_RELEASE_ASSERT(impl_->fr != nullptr, "Frame not available when accessed");
 478     return *impl_->fr;
 479 }
 480
 481 } // namespace gmx