src/gromacs/trajectoryanalysis/modules/pairdist.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2014,2015,2016,2018,2019, 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::PairDistance.
  38  *
  39  * \author Teemu Murtola <teemu.murtola@gmail.com>
  40  * \ingroup module_trajectoryanalysis
  41  */
  42 #include "gmxpre.h"
  43
  44 #include "pairdist.h"
  45
  46 #include <cmath>
  47
  48 #include <algorithm>
  49 #include <limits>
  50 #include <string>
  51 #include <vector>
  52
  53 #include "gromacs/analysisdata/analysisdata.h"
  54 #include "gromacs/analysisdata/modules/plot.h"
  55 #include "gromacs/options/basicoptions.h"
  56 #include "gromacs/options/filenameoption.h"
  57 #include "gromacs/options/ioptionscontainer.h"
  58 #include "gromacs/selection/nbsearch.h"
  59 #include "gromacs/selection/selection.h"
  60 #include "gromacs/selection/selectionoption.h"
  61 #include "gromacs/trajectory/trajectoryframe.h"
  62 #include "gromacs/trajectoryanalysis/analysissettings.h"
  63 #include "gromacs/trajectoryanalysis/topologyinformation.h"
  64 #include "gromacs/utility/arrayref.h"
  65 #include "gromacs/utility/exceptions.h"
  66 #include "gromacs/utility/stringutil.h"
  67
  68 namespace gmx
  69 {
  70
  71 namespace analysismodules
  72 {
  73
  74 namespace
  75 {
  76
  77 //! \addtogroup module_trajectoryanalysis
  78 //! \{
  79
  80 //! Enum value to store the selected value for `-type`.
  81 enum DistanceType
  82 {
  83     eDistanceType_Min,
  84     eDistanceType_Max
  85 };
  86
  87 //! Enum value to store the selected value for `-refgrouping`/`-selgrouping`.
  88 enum GroupType
  89 {
  90     eGroupType_All,
  91     eGroupType_Residue,
  92     eGroupType_Molecule,
  93     eGroupType_None
  94 };
  95
  96 //! Strings corresponding to DistanceType.
  97 const char* const c_distanceTypes[] = { "min", "max" };
  98 //! Strings corresponding to GroupType.
  99 const char* const c_groupTypes[] = { "all", "res", "mol", "none" };
 100
 101 /*! \brief
 102  * Implements `gmx pairdist` trajectory analysis module.
 103  */
 104 class PairDistance : public TrajectoryAnalysisModule
 105 {
 106 public:
 107     PairDistance();
 108
 109     void initOptions(IOptionsContainer* options, TrajectoryAnalysisSettings* settings) override;
 110     void initAnalysis(const TrajectoryAnalysisSettings& settings, const TopologyInformation& top) override;
 111
 112     TrajectoryAnalysisModuleDataPointer startFrames(const AnalysisDataParallelOptions& opt,
 113                                                     const SelectionCollection& selections) override;
 114     void                                analyzeFrame(int frnr, const t_trxframe& fr, t_pbc* pbc, TrajectoryAnalysisModuleData* pdata) override;
 115
 116     void finishAnalysis(int nframes) override;
 117     void writeOutput() override;
 118
 119 private:
 120     /*! \brief
 121      * Computed distances as a function of time.
 122      *
 123      * There is one data set for each selection in `sel_`.
 124      * Within each data set, there is one column for each distance to be
 125      * computed, as explained in the `-h` text.
 126      */
 127     AnalysisData distances_;
 128
 129     /*! \brief
 130      * Reference selection to compute distances to.
 131      *
 132      * mappedId() identifies the group (of type `refGroupType_`) into which
 133      * each position belogs.
 134      */
 135     Selection refSel_;
 136     /*! \brief
 137      * Selections to compute distances from.
 138      *
 139      * mappedId() identifies the group (of type `selGroupType_`) into which
 140      * each position belogs.
 141      */
 142     SelectionList sel_;
 143
 144     std::string fnDist_;
 145
 146     double       cutoff_;
 147     DistanceType distanceType_;
 148     GroupType    refGroupType_;
 149     GroupType    selGroupType_;
 150
 151     //! Number of groups in `refSel_`.
 152     int refGroupCount_;
 153     //! Maximum number of pairs of groups for one selection.
 154     int maxGroupCount_;
 155     //! Initial squared distance for distance accumulation.
 156     real initialDist2_;
 157     //! Cutoff squared for use in the actual calculation.
 158     real cutoff2_;
 159
 160     //! Neighborhood search object for the pair search.
 161     AnalysisNeighborhood nb_;
 162
 163     // Copy and assign disallowed by base.
 164 };
 165
 166 PairDistance::PairDistance() :
 167     cutoff_(0.0),
 168     distanceType_(eDistanceType_Min),
 169     refGroupType_(eGroupType_All),
 170     selGroupType_(eGroupType_All),
 171     refGroupCount_(0),
 172     maxGroupCount_(0),
 173     initialDist2_(0.0),
 174     cutoff2_(0.0)
 175 {
 176     registerAnalysisDataset(&distances_, "dist");
 177 }
 178
 179
 180 void PairDistance::initOptions(IOptionsContainer* options, TrajectoryAnalysisSettings* settings)
 181 {
 182     static const char* const desc[] = {
 183         "[THISMODULE] calculates pairwise distances between one reference",
 184         "selection (given with [TT]-ref[tt]) and one or more other selections",
 185         "(given with [TT]-sel[tt]).  It can calculate either the minimum",
 186         "distance (the default), or the maximum distance (with",
 187         "[TT]-type max[tt]).  Distances to each selection provided with",
 188         "[TT]-sel[tt] are computed independently.[PAR]",
 189         "By default, the global minimum/maximum distance is computed.",
 190         "To compute more distances (e.g., minimum distances to each residue",
 191         "in [TT]-ref[tt]), use [TT]-refgrouping[tt] and/or [TT]-selgrouping[tt]",
 192         "to specify how the positions within each selection should be",
 193         "grouped.[PAR]",
 194         "Computed distances are written to the file specified with [TT]-o[tt].",
 195         "If there are N groups in [TT]-ref[tt] and M groups in the first",
 196         "selection in [TT]-sel[tt], then the output contains N*M columns",
 197         "for the first selection. The columns contain distances like this:",
 198         "r1-s1, r2-s1, ..., r1-s2, r2-s2, ..., where rn is the n'th group",
 199         "in [TT]-ref[tt] and sn is the n'th group in the other selection.",
 200         "The distances for the second selection comes as separate columns",
 201         "after the first selection, and so on.  If some selections are",
 202         "dynamic, only the selected positions are used in the computation",
 203         "but the same number of columns is always written out.  If there",
 204         "are no positions contributing to some group pair, then the cutoff",
 205         "value is written (see below).[PAR]",
 206         "[TT]-cutoff[tt] sets a cutoff for the computed distances.",
 207         "If the result would contain a distance over the cutoff, the cutoff",
 208         "value is written to the output file instead. By default, no cutoff",
 209         "is used, but if you are not interested in values beyond a cutoff,",
 210         "or if you know that the minimum distance is smaller than a cutoff,",
 211         "you should set this option to allow the tool to use grid-based",
 212         "searching and be significantly faster.[PAR]",
 213         "If you want to compute distances between fixed pairs,",
 214         "[gmx-distance] may be a more suitable tool."
 215     };
 216
 217     settings->setHelpText(desc);
 218
 219     options->addOption(FileNameOption("o")
 220                                .filetype(eftPlot)
 221                                .outputFile()
 222                                .required()
 223                                .store(&fnDist_)
 224                                .defaultBasename("dist")
 225                                .description("Distances as function of time"));
 226
 227     options->addOption(
 228             DoubleOption("cutoff").store(&cutoff_).description("Maximum distance to consider"));
 229     options->addOption(EnumOption<DistanceType>("type")
 230                                .store(&distanceType_)
 231                                .enumValue(c_distanceTypes)
 232                                .description("Type of distances to calculate"));
 233     options->addOption(
 234             EnumOption<GroupType>("refgrouping")
 235                     .store(&refGroupType_)
 236                     .enumValue(c_groupTypes)
 237                     .description("Grouping of -ref positions to compute the min/max over"));
 238     options->addOption(
 239             EnumOption<GroupType>("selgrouping")
 240                     .store(&selGroupType_)
 241                     .enumValue(c_groupTypes)
 242                     .description("Grouping of -sel positions to compute the min/max over"));
 243
 244     options->addOption(SelectionOption("ref").store(&refSel_).required().description(
 245             "Reference positions to calculate distances from"));
 246     options->addOption(SelectionOption("sel").storeVector(&sel_).required().multiValue().description(
 247             "Positions to calculate distances for"));
 248 }
 249
 250 //! Helper function to initialize the grouping for a selection.
 251 int initSelectionGroups(Selection* sel, const gmx_mtop_t* top, int type)
 252 {
 253     e_index_t indexType = INDEX_UNKNOWN;
 254     switch (type)
 255     {
 256         case eGroupType_All: indexType = INDEX_ALL; break;
 257         case eGroupType_Residue: indexType = INDEX_RES; break;
 258         case eGroupType_Molecule: indexType = INDEX_MOL; break;
 259         case eGroupType_None: indexType = INDEX_ATOM; break;
 260     }
 261     return sel->initOriginalIdsToGroup(top, indexType);
 262 }
 263
 264
 265 void PairDistance::initAnalysis(const TrajectoryAnalysisSettings& settings, const TopologyInformation& top)
 266 {
 267     refGroupCount_ = initSelectionGroups(&refSel_, top.mtop(), refGroupType_);
 268
 269     maxGroupCount_ = 0;
 270     distances_.setDataSetCount(sel_.size());
 271     for (size_t i = 0; i < sel_.size(); ++i)
 272     {
 273         const int selGroupCount = initSelectionGroups(&sel_[i], top.mtop(), selGroupType_);
 274         const int columnCount   = refGroupCount_ * selGroupCount;
 275         maxGroupCount_          = std::max(maxGroupCount_, columnCount);
 276         distances_.setColumnCount(i, columnCount);
 277     }
 278
 279     if (!fnDist_.empty())
 280     {
 281         AnalysisDataPlotModulePointer plotm(new AnalysisDataPlotModule(settings.plotSettings()));
 282         plotm->setFileName(fnDist_);
 283         if (distanceType_ == eDistanceType_Max)
 284         {
 285             plotm->setTitle("Maximum distance");
 286         }
 287         else
 288         {
 289             plotm->setTitle("Minimum distance");
 290         }
 291         // TODO: Figure out and add a descriptive subtitle and/or a longer
 292         // title and/or better legends based on the grouping and the reference
 293         // selection.
 294         plotm->setXAxisIsTime();
 295         plotm->setYLabel("Distance (nm)");
 296         for (size_t g = 0; g < sel_.size(); ++g)
 297         {
 298             plotm->appendLegend(sel_[g].name());
 299         }
 300         distances_.addModule(plotm);
 301     }
 302
 303     nb_.setCutoff(cutoff_);
 304     if (cutoff_ <= 0.0)
 305     {
 306         cutoff_       = 0.0;
 307         initialDist2_ = std::numeric_limits<real>::max();
 308     }
 309     else
 310     {
 311         initialDist2_ = cutoff_ * cutoff_;
 312     }
 313     if (distanceType_ == eDistanceType_Max)
 314     {
 315         initialDist2_ = 0.0;
 316     }
 317     cutoff2_ = cutoff_ * cutoff_;
 318 }
 319
 320 /*! \brief
 321  * Temporary memory for use within a single-frame calculation.
 322  */
 323 class PairDistanceModuleData : public TrajectoryAnalysisModuleData
 324 {
 325 public:
 326     /*! \brief
 327      * Reserves memory for the frame-local data.
 328      */
 329     PairDistanceModuleData(TrajectoryAnalysisModule*          module,
 330                            const AnalysisDataParallelOptions& opt,
 331                            const SelectionCollection&         selections,
 332                            int                                refGroupCount,
 333                            const Selection&                   refSel,
 334                            int                                maxGroupCount) :
 335         TrajectoryAnalysisModuleData(module, opt, selections)
 336     {
 337         distArray_.resize(maxGroupCount);
 338         countArray_.resize(maxGroupCount);
 339         refCountArray_.resize(refGroupCount);
 340         if (!refSel.isDynamic())
 341         {
 342             initRefCountArray(refSel);
 343         }
 344     }
 345
 346     void finish() override { finishDataHandles(); }
 347
 348     /*! \brief
 349      * Computes the number of positions in each group in \p refSel
 350      * and stores them into `refCountArray_`.
 351      */
 352     void initRefCountArray(const Selection& refSel)
 353     {
 354         std::fill(refCountArray_.begin(), refCountArray_.end(), 0);
 355         int refPos = 0;
 356         while (refPos < refSel.posCount())
 357         {
 358             const int refIndex = refSel.position(refPos).mappedId();
 359             const int startPos = refPos;
 360             ++refPos;
 361             while (refPos < refSel.posCount() && refSel.position(refPos).mappedId() == refIndex)
 362             {
 363                 ++refPos;
 364             }
 365             refCountArray_[refIndex] = refPos - startPos;
 366         }
 367     }
 368
 369     /*! \brief
 370      * Squared distance between each group
 371      *
 372      * One entry for each group pair for the current selection.
 373      * Enough memory is allocated to fit the largest calculation selection.
 374      * This is needed to support neighborhood searching, which may not
 375      * return the pairs in order: for each group pair, we need to search
 376      * through all the position pairs and update this array to find the
 377      * minimum/maximum distance between them.
 378      */
 379     std::vector<real> distArray_;
 380     /*! \brief
 381      * Number of pairs within the cutoff that have contributed to the value
 382      * in `distArray_`.
 383      *
 384      * This is needed to identify whether there were any pairs inside the
 385      * cutoff and whether there were additional pairs outside the cutoff
 386      * that were not covered by the neihborhood search.
 387      */
 388     std::vector<int> countArray_;
 389     /*! \brief
 390      * Number of positions within each reference group.
 391      *
 392      * This is used to more efficiently compute the total number of pairs
 393      * (for comparison with `countArray_`), as otherwise these numbers
 394      * would need to be recomputed for each selection.
 395      */
 396     std::vector<int> refCountArray_;
 397 };
 398
 399 TrajectoryAnalysisModuleDataPointer PairDistance::startFrames(const AnalysisDataParallelOptions& opt,
 400                                                               const SelectionCollection& selections)
 401 {
 402     return TrajectoryAnalysisModuleDataPointer(new PairDistanceModuleData(
 403             this, opt, selections, refGroupCount_, refSel_, maxGroupCount_));
 404 }
 405
 406 void PairDistance::analyzeFrame(int frnr, const t_trxframe& fr, t_pbc* pbc, TrajectoryAnalysisModuleData* pdata)
 407 {
 408     AnalysisDataHandle      dh         = pdata->dataHandle(distances_);
 409     const Selection&        refSel     = pdata->parallelSelection(refSel_);
 410     const SelectionList&    sel        = pdata->parallelSelections(sel_);
 411     PairDistanceModuleData& frameData  = *static_cast<PairDistanceModuleData*>(pdata);
 412     std::vector<real>&      distArray  = frameData.distArray_;
 413     std::vector<int>&       countArray = frameData.countArray_;
 414
 415     if (cutoff_ > 0.0 && refSel.isDynamic())
 416     {
 417         // Count the number of reference positions in each group, so that
 418         // this does not need to be computed again for each selection.
 419         // This is needed only if it is possible that the neighborhood search
 420         // does not cover all the pairs, hence the cutoff > 0.0 check.
 421         // If refSel is static, then the array contents are static as well,
 422         // and it has been initialized in the constructor of the data object.
 423         frameData.initRefCountArray(refSel);
 424     }
 425     const std::vector<int>& refCountArray = frameData.refCountArray_;
 426
 427     AnalysisNeighborhoodSearch nbsearch = nb_.initSearch(pbc, refSel);
 428     dh.startFrame(frnr, fr.time);
 429     for (size_t g = 0; g < sel.size(); ++g)
 430     {
 431         const int columnCount = distances_.columnCount(g);
 432         std::fill(distArray.begin(), distArray.begin() + columnCount, initialDist2_);
 433         std::fill(countArray.begin(), countArray.begin() + columnCount, 0);
 434
 435         // Accumulate the number of position pairs within the cutoff and the
 436         // min/max distance for each group pair.
 437         AnalysisNeighborhoodPairSearch pairSearch = nbsearch.startPairSearch(sel[g]);
 438         AnalysisNeighborhoodPair       pair;
 439         while (pairSearch.findNextPair(&pair))
 440         {
 441             const SelectionPosition& refPos   = refSel.position(pair.refIndex());
 442             const SelectionPosition& selPos   = sel[g].position(pair.testIndex());
 443             const int                refIndex = refPos.mappedId();
 444             const int                selIndex = selPos.mappedId();
 445             const int                index    = selIndex * refGroupCount_ + refIndex;
 446             const real               r2       = pair.distance2();
 447             if (distanceType_ == eDistanceType_Min)
 448             {
 449                 if (distArray[index] > r2)
 450                 {
 451                     distArray[index] = r2;
 452                 }
 453             }
 454             else
 455             {
 456                 if (distArray[index] < r2)
 457                 {
 458                     distArray[index] = r2;
 459                 }
 460             }
 461             ++countArray[index];
 462         }
 463
 464         // If it is possible that positions outside the cutoff (or lack of
 465         // them) affects the result, then we need to check whether there were
 466         // any.  This is necessary for two cases:
 467         //  - With max distances, if there are pairs outside the cutoff, then
 468         //    the computed distance should be equal to the cutoff instead of
 469         //    the largest distance that was found above.
 470         //  - With either distance type, if all pairs are outside the cutoff,
 471         //    then countArray must be updated so that the presence flag
 472         //    in the output data reflects the dynamic selection status, not
 473         //    whether something was inside the cutoff or not.
 474         if (cutoff_ > 0.0)
 475         {
 476             int selPos = 0;
 477             // Loop over groups in this selection (at start, selPos is always
 478             // the first position in the next group).
 479             while (selPos < sel[g].posCount())
 480             {
 481                 // Count the number of positions in this group.
 482                 const int selIndex = sel[g].position(selPos).mappedId();
 483                 const int startPos = selPos;
 484                 ++selPos;
 485                 while (selPos < sel[g].posCount() && sel[g].position(selPos).mappedId() == selIndex)
 486                 {
 487                     ++selPos;
 488                 }
 489                 const int count = selPos - startPos;
 490                 // Check all group pairs that contain this group.
 491                 for (int i = 0; i < refGroupCount_; ++i)
 492                 {
 493                     const int index      = selIndex * refGroupCount_ + i;
 494                     const int totalCount = refCountArray[i] * count;
 495                     // If there were positions outside the cutoff,
 496                     // update the distance if necessary and the count.
 497                     if (countArray[index] < totalCount)
 498                     {
 499                         if (distanceType_ == eDistanceType_Max)
 500                         {
 501                             distArray[index] = cutoff2_;
 502                         }
 503                         countArray[index] = totalCount;
 504                     }
 505                 }
 506             }
 507         }
 508
 509         // Write the computed distances to the output data.
 510         dh.selectDataSet(g);
 511         for (int i = 0; i < columnCount; ++i)
 512         {
 513             if (countArray[i] > 0)
 514             {
 515                 dh.setPoint(i, std::sqrt(distArray[i]));
 516             }
 517             else
 518             {
 519                 // If there are no contributing positions, write out the cutoff
 520                 // value.
 521                 dh.setPoint(i, cutoff_, false);
 522             }
 523         }
 524     }
 525     dh.finishFrame();
 526 }
 527
 528 void PairDistance::finishAnalysis(int /*nframes*/) {}
 529
 530 void PairDistance::writeOutput() {}
 531
 532 //! \}
 533
 534 } // namespace
 535
 536 const char PairDistanceInfo::name[] = "pairdist";
 537 const char PairDistanceInfo::shortDescription[] =
 538         "Calculate pairwise distances between groups of positions";
 539
 540 TrajectoryAnalysisModulePointer PairDistanceInfo::create()
 541 {
 542     return TrajectoryAnalysisModulePointer(new PairDistance);
 543 }
 544
 545 } // namespace analysismodules
 546
 547 } // namespace gmx