src/gromacs/trajectoryanalysis/modules/pairdist.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2014,2015, 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/fileio/trx.h"
  56 #include "gromacs/options/basicoptions.h"
  57 #include "gromacs/options/filenameoption.h"
  58 #include "gromacs/options/ioptionscontainer.h"
  59 #include "gromacs/selection/nbsearch.h"
  60 #include "gromacs/selection/selection.h"
  61 #include "gromacs/selection/selectionoption.h"
  62 #include "gromacs/trajectoryanalysis/analysissettings.h"
  63 #include "gromacs/utility/arrayref.h"
  64 #include "gromacs/utility/exceptions.h"
  65 #include "gromacs/utility/stringutil.h"
  66
  67 namespace gmx
  68 {
  69
  70 namespace analysismodules
  71 {
  72
  73 namespace
  74 {
  75
  76 //! \addtogroup module_trajectoryanalysis
  77 //! \{
  78
  79 //! Enum value to store the selected value for `-type`.
  80 enum DistanceType
  81 {
  82     eDistanceType_Min,
  83     eDistanceType_Max
  84 };
  85
  86 //! Enum value to store the selected value for `-refgrouping`/`-selgrouping`.
  87 enum GroupType
  88 {
  89     eGroupType_All,
  90     eGroupType_Residue,
  91     eGroupType_Molecule,
  92     eGroupType_None
  93 };
  94
  95 //! Strings corresponding to DistanceType.
  96 const char *const           c_distanceTypes[] = { "min", "max" };
  97 //! Strings corresponding to GroupType.
  98 const char *const           c_groupTypes[]    = { "all", "res", "mol", "none" };
  99
 100 /*! \brief
 101  * Implements `gmx pairdist` trajectory analysis module.
 102  */
 103 class PairDistance : public TrajectoryAnalysisModule
 104 {
 105     public:
 106         PairDistance();
 107
 108         virtual void initOptions(IOptionsContainer          *options,
 109                                  TrajectoryAnalysisSettings *settings);
 110         virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
 111                                   const TopologyInformation        &top);
 112
 113         virtual TrajectoryAnalysisModuleDataPointer startFrames(
 114             const AnalysisDataParallelOptions &opt,
 115             const SelectionCollection         &selections);
 116         virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 117                                   TrajectoryAnalysisModuleData *pdata);
 118
 119         virtual void finishAnalysis(int nframes);
 120         virtual void writeOutput();
 121
 122     private:
 123         /*! \brief
 124          * Computed distances as a function of time.
 125          *
 126          * There is one data set for each selection in `sel_`.
 127          * Within each data set, there is one column for each distance to be
 128          * computed, as explained in the `-h` text.
 129          */
 130         AnalysisData            distances_;
 131
 132         /*! \brief
 133          * Reference selection to compute distances to.
 134          *
 135          * mappedId() identifies the group (of type `refGroupType_`) into which
 136          * each position belogs.
 137          */
 138         Selection               refSel_;
 139         /*! \brief
 140          * Selections to compute distances from.
 141          *
 142          * mappedId() identifies the group (of type `selGroupType_`) into which
 143          * each position belogs.
 144          */
 145         SelectionList           sel_;
 146
 147         std::string             fnDist_;
 148
 149         double                  cutoff_;
 150         int                     distanceType_;
 151         int                     refGroupType_;
 152         int                     selGroupType_;
 153
 154         //! Number of groups in `refSel_`.
 155         int                     refGroupCount_;
 156         //! Maximum number of pairs of groups for one selection.
 157         int                     maxGroupCount_;
 158         //! Initial squared distance for distance accumulation.
 159         real                    initialDist2_;
 160         //! Cutoff squared for use in the actual calculation.
 161         real                    cutoff2_;
 162
 163         //! Neighborhood search object for the pair search.
 164         AnalysisNeighborhood    nb_;
 165
 166         // Copy and assign disallowed by base.
 167 };
 168
 169 PairDistance::PairDistance()
 170     : TrajectoryAnalysisModule(PairDistanceInfo::name, PairDistanceInfo::shortDescription),
 171       cutoff_(0.0), distanceType_(eDistanceType_Min),
 172       refGroupType_(eGroupType_All), selGroupType_(eGroupType_All),
 173       refGroupCount_(0), maxGroupCount_(0), initialDist2_(0.0), cutoff2_(0.0)
 174 {
 175     registerAnalysisDataset(&distances_, "dist");
 176 }
 177
 178
 179 void
 180 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").filetype(eftPlot).outputFile().required()
 220                            .store(&fnDist_).defaultBasename("dist")
 221                            .description("Distances as function of time"));
 222
 223     options->addOption(DoubleOption("cutoff").store(&cutoff_)
 224                            .description("Maximum distance to consider"));
 225     options->addOption(StringOption("type").storeEnumIndex(&distanceType_)
 226                            .defaultEnumIndex(0).enumValue(c_distanceTypes)
 227                            .description("Type of distances to calculate"));
 228     options->addOption(StringOption("refgrouping").storeEnumIndex(&refGroupType_)
 229                            .defaultEnumIndex(0).enumValue(c_groupTypes)
 230                            .description("Grouping of -ref positions to compute the min/max over"));
 231     options->addOption(StringOption("selgrouping").storeEnumIndex(&selGroupType_)
 232                            .defaultEnumIndex(0).enumValue(c_groupTypes)
 233                            .description("Grouping of -sel positions to compute the min/max over"));
 234
 235     options->addOption(SelectionOption("ref").store(&refSel_).required()
 236                            .description("Reference positions to calculate distances from"));
 237     options->addOption(SelectionOption("sel").storeVector(&sel_).required().multiValue()
 238                            .description("Positions to calculate distances for"));
 239 }
 240
 241 //! Helper function to initialize the grouping for a selection.
 242 int initSelectionGroups(Selection *sel, t_topology *top, int type)
 243 {
 244     e_index_t indexType = INDEX_UNKNOWN;
 245     switch (type)
 246     {
 247         case eGroupType_All:      indexType = INDEX_ALL; break;
 248         case eGroupType_Residue:  indexType = INDEX_RES; break;
 249         case eGroupType_Molecule: indexType = INDEX_MOL; break;
 250         case eGroupType_None:     indexType = INDEX_ATOM; break;
 251     }
 252     return sel->initOriginalIdsToGroup(top, indexType);
 253 }
 254
 255
 256 void
 257 PairDistance::initAnalysis(const TrajectoryAnalysisSettings &settings,
 258                            const TopologyInformation        &top)
 259 {
 260     refGroupCount_ = initSelectionGroups(&refSel_, top.topology(), refGroupType_);
 261
 262     maxGroupCount_ = 0;
 263     distances_.setDataSetCount(sel_.size());
 264     for (size_t i = 0; i < sel_.size(); ++i)
 265     {
 266         const int selGroupCount
 267             = initSelectionGroups(&sel_[i], top.topology(), selGroupType_);
 268         const int columnCount = refGroupCount_ * selGroupCount;
 269         maxGroupCount_ = std::max(maxGroupCount_, columnCount);
 270         distances_.setColumnCount(i, columnCount);
 271     }
 272
 273     if (!fnDist_.empty())
 274     {
 275         AnalysisDataPlotModulePointer plotm(
 276                 new AnalysisDataPlotModule(settings.plotSettings()));
 277         plotm->setFileName(fnDist_);
 278         if (distanceType_ == eDistanceType_Max)
 279         {
 280             plotm->setTitle("Maximum distance");
 281         }
 282         else
 283         {
 284             plotm->setTitle("Minimum distance");
 285         }
 286         // TODO: Figure out and add a descriptive subtitle and/or a longer
 287         // title and/or better legends based on the grouping and the reference
 288         // selection.
 289         plotm->setXAxisIsTime();
 290         plotm->setYLabel("Distance (nm)");
 291         for (size_t g = 0; g < sel_.size(); ++g)
 292         {
 293             plotm->appendLegend(sel_[g].name());
 294         }
 295         distances_.addModule(plotm);
 296     }
 297
 298     nb_.setCutoff(cutoff_);
 299     if (cutoff_ <= 0.0)
 300     {
 301         cutoff_       = 0.0;
 302         initialDist2_ = std::numeric_limits<real>::max();
 303     }
 304     else
 305     {
 306         initialDist2_ = cutoff_ * cutoff_;
 307     }
 308     if (distanceType_ == eDistanceType_Max)
 309     {
 310         initialDist2_ = 0.0;
 311     }
 312     cutoff2_ = cutoff_ * cutoff_;
 313 }
 314
 315 /*! \brief
 316  * Temporary memory for use within a single-frame calculation.
 317  */
 318 class PairDistanceModuleData : public TrajectoryAnalysisModuleData
 319 {
 320     public:
 321         /*! \brief
 322          * Reserves memory for the frame-local data.
 323          */
 324         PairDistanceModuleData(TrajectoryAnalysisModule          *module,
 325                                const AnalysisDataParallelOptions &opt,
 326                                const SelectionCollection         &selections,
 327                                int                                refGroupCount,
 328                                const Selection                   &refSel,
 329                                int                                maxGroupCount)
 330             : TrajectoryAnalysisModuleData(module, opt, selections)
 331         {
 332             distArray_.resize(maxGroupCount);
 333             countArray_.resize(maxGroupCount);
 334             refCountArray_.resize(refGroupCount);
 335             if (!refSel.isDynamic())
 336             {
 337                 initRefCountArray(refSel);
 338             }
 339         }
 340
 341         virtual void finish() { finishDataHandles(); }
 342
 343         /*! \brief
 344          * Computes the number of positions in each group in \p refSel
 345          * and stores them into `refCountArray_`.
 346          */
 347         void initRefCountArray(const Selection &refSel)
 348         {
 349             std::fill(refCountArray_.begin(), refCountArray_.end(), 0);
 350             int refPos = 0;
 351             while (refPos < refSel.posCount())
 352             {
 353                 const int refIndex = refSel.position(refPos).mappedId();
 354                 const int startPos = refPos;
 355                 ++refPos;
 356                 while (refPos < refSel.posCount()
 357                        && refSel.position(refPos).mappedId() == refIndex)
 358                 {
 359                     ++refPos;
 360                 }
 361                 refCountArray_[refIndex] = refPos - startPos;
 362             }
 363         }
 364
 365         /*! \brief
 366          * Squared distance between each group
 367          *
 368          * One entry for each group pair for the current selection.
 369          * Enough memory is allocated to fit the largest calculation selection.
 370          * This is needed to support neighborhood searching, which may not
 371          * return the pairs in order: for each group pair, we need to search
 372          * through all the position pairs and update this array to find the
 373          * minimum/maximum distance between them.
 374          */
 375         std::vector<real> distArray_;
 376         /*! \brief
 377          * Number of pairs within the cutoff that have contributed to the value
 378          * in `distArray_`.
 379          *
 380          * This is needed to identify whether there were any pairs inside the
 381          * cutoff and whether there were additional pairs outside the cutoff
 382          * that were not covered by the neihborhood search.
 383          */
 384         std::vector<int>  countArray_;
 385         /*! \brief
 386          * Number of positions within each reference group.
 387          *
 388          * This is used to more efficiently compute the total number of pairs
 389          * (for comparison with `countArray_`), as otherwise these numbers
 390          * would need to be recomputed for each selection.
 391          */
 392         std::vector<int>  refCountArray_;
 393 };
 394
 395 TrajectoryAnalysisModuleDataPointer PairDistance::startFrames(
 396         const AnalysisDataParallelOptions &opt,
 397         const SelectionCollection         &selections)
 398 {
 399     return TrajectoryAnalysisModuleDataPointer(
 400             new PairDistanceModuleData(this, opt, selections, refGroupCount_,
 401                                        refSel_, maxGroupCount_));
 402 }
 403
 404 void
 405 PairDistance::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 406                            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()
 486                        && sel[g].position(selPos).mappedId() == selIndex)
 487                 {
 488                     ++selPos;
 489                 }
 490                 const int count = selPos - startPos;
 491                 // Check all group pairs that contain this group.
 492                 for (int i = 0; i < refGroupCount_; ++i)
 493                 {
 494                     const int index      = selIndex * refGroupCount_ + i;
 495                     const int totalCount = refCountArray[i] * count;
 496                     // If there were positions outside the cutoff,
 497                     // update the distance if necessary and the count.
 498                     if (countArray[index] < totalCount)
 499                     {
 500                         if (distanceType_ == eDistanceType_Max)
 501                         {
 502                             distArray[index] = cutoff2_;
 503                         }
 504                         countArray[index] = totalCount;
 505                     }
 506                 }
 507             }
 508         }
 509
 510         // Write the computed distances to the output data.
 511         dh.selectDataSet(g);
 512         for (int i = 0; i < columnCount; ++i)
 513         {
 514             if (countArray[i] > 0)
 515             {
 516                 dh.setPoint(i, std::sqrt(distArray[i]));
 517             }
 518             else
 519             {
 520                 // If there are no contributing positions, write out the cutoff
 521                 // value.
 522                 dh.setPoint(i, cutoff_, false);
 523             }
 524         }
 525     }
 526     dh.finishFrame();
 527 }
 528
 529 void
 530 PairDistance::finishAnalysis(int /*nframes*/)
 531 {
 532 }
 533
 534 void
 535 PairDistance::writeOutput()
 536 {
 537 }
 538
 539 //! \}
 540
 541 }       // namespace
 542
 543 const char PairDistanceInfo::name[]             = "pairdist";
 544 const char PairDistanceInfo::shortDescription[] =
 545     "Calculate pairwise distances between groups of positions";
 546
 547 TrajectoryAnalysisModulePointer PairDistanceInfo::create()
 548 {
 549     return TrajectoryAnalysisModulePointer(new PairDistance);
 550 }
 551
 552 } // namespace analysismodules
 553
 554 } // namespace gmx