src/gromacs/selection/nbsearch.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2009,2010,2011,2012,2013,2014, 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 neighborhood searching for analysis (from nbsearch.h).
  38  *
  39  * \todo
  40  * The grid implementation could still be optimized in several different ways:
  41  *   - Triclinic grid cells are not the most efficient shape, but make PBC
  42  *     handling easier.
  43  *   - Precalculating the required PBC shift for a pair of cells outside the
  44  *     inner loop. After this is done, it should be quite straightforward to
  45  *     move to rectangular cells.
  46  *   - Pruning grid cells from the search list if they are completely outside
  47  *     the sphere that is being considered.
  48  *   - A better heuristic could be added for falling back to simple loops for a
  49  *     small number of reference particles.
  50  *   - A better heuristic for selecting the grid size.
  51  *   - A multi-level grid implementation could be used to be able to use small
  52  *     grids for short cutoffs with very inhomogeneous particle distributions
  53  *     without a memory cost.
  54  *
  55  * \author Teemu Murtola <teemu.murtola@gmail.com>
  56  * \ingroup module_selection
  57  */
  58 #include "gmxpre.h"
  59
  60 #include "gromacs/selection/nbsearch.h"
  61
  62 #include <cmath>
  63 #include <cstring>
  64
  65 #include <algorithm>
  66 #include <vector>
  67
  68 #include "thread_mpi/mutex.h"
  69
  70 #include "gromacs/legacyheaders/names.h"
  71
  72 #include "gromacs/math/vec.h"
  73 #include "gromacs/pbcutil/pbc.h"
  74 #include "gromacs/selection/position.h"
  75 #include "gromacs/topology/block.h"
  76 #include "gromacs/utility/arrayref.h"
  77 #include "gromacs/utility/exceptions.h"
  78 #include "gromacs/utility/gmxassert.h"
  79 #include "gromacs/utility/smalloc.h"
  80 #include "gromacs/utility/stringutil.h"
  81
  82 namespace gmx
  83 {
  84
  85 namespace internal
  86 {
  87
  88 /********************************************************************
  89  * Implementation class declarations
  90  */
  91
  92 class AnalysisNeighborhoodSearchImpl
  93 {
  94     public:
  95         typedef AnalysisNeighborhoodPairSearch::ImplPointer
  96             PairSearchImplPointer;
  97         typedef std::vector<PairSearchImplPointer> PairSearchList;
  98         typedef std::vector<std::vector<int> > CellList;
  99
 100         explicit AnalysisNeighborhoodSearchImpl(real cutoff);
 101         ~AnalysisNeighborhoodSearchImpl();
 102
 103         /*! \brief
 104          * Initializes the search with a given box and reference positions.
 105          *
 106          * \param[in]     mode      Search mode to use.
 107          * \param[in]     bXY       Whether to use 2D searching.
 108          * \param[in]     excls     Exclusions.
 109          * \param[in]     pbc       PBC information.
 110          * \param[in]     positions Set of reference positions.
 111          */
 112         void init(AnalysisNeighborhood::SearchMode     mode,
 113                   bool                                 bXY,
 114                   const t_blocka                      *excls,
 115                   const t_pbc                         *pbc,
 116                   const AnalysisNeighborhoodPositions &positions);
 117         PairSearchImplPointer getPairSearch();
 118
 119         real cutoffSquared() const { return cutoff2_; }
 120         bool usesGridSearch() const { return bGrid_; }
 121
 122     private:
 123         //! Calculates offsets to neighboring grid cells that should be considered.
 124         void initGridCellNeighborList();
 125         /*! \brief
 126          * Determines a suitable grid size and sets up the cells.
 127          *
 128          * \param[in]     pbc  Information about the box.
 129          * \returns  false if grid search is not suitable.
 130          */
 131         bool initGridCells(const t_pbc &pbc);
 132         /*! \brief
 133          * Sets ua a search grid for a given box.
 134          *
 135          * \param[in]     pbc  Information about the box.
 136          * \returns  false if grid search is not suitable.
 137          */
 138         bool initGrid(const t_pbc &pbc);
 139         /*! \brief
 140          * Maps a point into a grid cell.
 141          *
 142          * \param[in]  x    Point to map.
 143          * \param[out] cell Indices of the grid cell in which \p x lies.
 144          * \param[out] xout Coordinates to use
 145          *     (will be within the triclinic unit cell).
 146          */
 147         void mapPointToGridCell(const rvec x, ivec cell, rvec xout) const;
 148         /*! \brief
 149          * Calculates linear index of a grid cell.
 150          *
 151          * \param[in]  cell Cell indices.
 152          * \returns    Linear index of \p cell.
 153          */
 154         int getGridCellIndex(const ivec cell) const;
 155         /*! \brief
 156          * Adds an index into a grid cell.
 157          *
 158          * \param[in]     cell Cell into which \p i should be added.
 159          * \param[in]     i    Index to add.
 160          */
 161         void addToGridCell(const ivec cell, int i);
 162
 163         //! Whether to try grid searching.
 164         bool                    bTryGrid_;
 165         //! The cutoff.
 166         real                    cutoff_;
 167         //! The cutoff squared.
 168         real                    cutoff2_;
 169         //! Whether to do searching in XY plane only.
 170         bool                    bXY_;
 171
 172         //! Number of reference points for the current frame.
 173         int                     nref_;
 174         //! Reference point positions.
 175         const rvec             *xref_;
 176         //! Reference position exclusion IDs.
 177         const int              *refExclusionIds_;
 178         //! Exclusions.
 179         const t_blocka         *excls_;
 180         //! PBC data.
 181         t_pbc                   pbc_;
 182
 183         //! Whether grid searching is actually used for the current positions.
 184         bool                    bGrid_;
 185         //! Array allocated for storing in-unit-cell reference positions.
 186         rvec                   *xref_alloc_;
 187         //! Allocation count for xref_alloc.
 188         int                     xref_nalloc_;
 189         //! false if the box is rectangular.
 190         bool                    bTric_;
 191         //! Box vectors of a single grid cell.
 192         matrix                  cellbox_;
 193         //! The reciprocal cell vectors as columns; the inverse of \p cellbox.
 194         matrix                  recipcell_;
 195         //! Number of cells along each dimension.
 196         ivec                    ncelldim_;
 197         //! Data structure to hold the grid cell contents.
 198         CellList                cells_;
 199         //! Number of neighboring cells to consider.
 200         int                     ngridnb_;
 201         //! Offsets of the neighboring cells to consider.
 202         ivec                   *gnboffs_;
 203         //! Allocation count for \p gnboffs.
 204         int                     gnboffs_nalloc_;
 205
 206         tMPI::mutex             createPairSearchMutex_;
 207         PairSearchList          pairSearchList_;
 208
 209         friend class AnalysisNeighborhoodPairSearchImpl;
 210
 211         GMX_DISALLOW_COPY_AND_ASSIGN(AnalysisNeighborhoodSearchImpl);
 212 };
 213
 214 class AnalysisNeighborhoodPairSearchImpl
 215 {
 216     public:
 217         explicit AnalysisNeighborhoodPairSearchImpl(const AnalysisNeighborhoodSearchImpl &search)
 218             : search_(search)
 219         {
 220             testExclusionIds_ = NULL;
 221             nexcl_            = 0;
 222             excl_             = NULL;
 223             clear_rvec(xtest_);
 224             clear_ivec(testcell_);
 225             reset(-1);
 226         }
 227
 228         //! Initializes a search to find reference positions neighboring \p x.
 229         void startSearch(const AnalysisNeighborhoodPositions &positions);
 230         //! Searches for the next neighbor.
 231         template <class Action>
 232         bool searchNext(Action action);
 233         //! Initializes a pair representing the pair found by searchNext().
 234         void initFoundPair(AnalysisNeighborhoodPair *pair) const;
 235         //! Advances to the next test position, skipping any remaining pairs.
 236         void nextTestPosition();
 237
 238     private:
 239         //! Clears the loop indices.
 240         void reset(int testIndex);
 241         //! Checks whether a reference positiong should be excluded.
 242         bool isExcluded(int j);
 243
 244         //! Parent search object.
 245         const AnalysisNeighborhoodSearchImpl   &search_;
 246         //! Reference to the test positions.
 247         ConstArrayRef<rvec>                     testPositions_;
 248         //! Reference to the test exclusion indices.
 249         const int                              *testExclusionIds_;
 250         //! Number of excluded reference positions for current test particle.
 251         int                                     nexcl_;
 252         //! Exclusions for current test particle.
 253         const int                              *excl_;
 254         //! Index of the currently active test position in \p testPositions_.
 255         int                                     testIndex_;
 256         //! Stores test position during a pair loop.
 257         rvec                                    xtest_;
 258         //! Stores the previous returned position during a pair loop.
 259         int                                     previ_;
 260         //! Stores the pair distance corresponding to previ_;
 261         real                                    prevr2_;
 262         //! Stores the current exclusion index during loops.
 263         int                                     exclind_;
 264         //! Stores the test particle cell index during loops.
 265         ivec                                    testcell_;
 266         //! Stores the current cell neighbor index during pair loops.
 267         int                                     prevnbi_;
 268         //! Stores the index within the current cell during pair loops.
 269         int                                     prevcai_;
 270
 271         GMX_DISALLOW_COPY_AND_ASSIGN(AnalysisNeighborhoodPairSearchImpl);
 272 };
 273
 274 /********************************************************************
 275  * AnalysisNeighborhoodSearchImpl
 276  */
 277
 278 AnalysisNeighborhoodSearchImpl::AnalysisNeighborhoodSearchImpl(real cutoff)
 279 {
 280     bTryGrid_       = true;
 281     cutoff_         = cutoff;
 282     if (cutoff_ <= 0)
 283     {
 284         cutoff_     = GMX_REAL_MAX;
 285         bTryGrid_   = false;
 286     }
 287     cutoff2_        = sqr(cutoff_);
 288     bXY_            = false;
 289
 290     nref_            = 0;
 291     xref_            = NULL;
 292     refExclusionIds_ = NULL;
 293     std::memset(&pbc_, 0, sizeof(pbc_));
 294
 295     bGrid_          = false;
 296
 297     xref_alloc_     = NULL;
 298     xref_nalloc_    = 0;
 299     bTric_          = false;
 300     clear_mat(cellbox_);
 301     clear_mat(recipcell_);
 302     clear_ivec(ncelldim_);
 303
 304     ngridnb_        = 0;
 305     gnboffs_        = NULL;
 306     gnboffs_nalloc_ = 0;
 307 }
 308
 309 AnalysisNeighborhoodSearchImpl::~AnalysisNeighborhoodSearchImpl()
 310 {
 311     PairSearchList::const_iterator i;
 312     for (i = pairSearchList_.begin(); i != pairSearchList_.end(); ++i)
 313     {
 314         GMX_RELEASE_ASSERT(i->unique(),
 315                            "Dangling AnalysisNeighborhoodPairSearch reference");
 316     }
 317     sfree(xref_alloc_);
 318     sfree(gnboffs_);
 319 }
 320
 321 AnalysisNeighborhoodSearchImpl::PairSearchImplPointer
 322 AnalysisNeighborhoodSearchImpl::getPairSearch()
 323 {
 324     tMPI::lock_guard<tMPI::mutex> lock(createPairSearchMutex_);
 325     // TODO: Consider whether this needs to/can be faster, e.g., by keeping a
 326     // separate pool of unused search objects.
 327     PairSearchList::const_iterator i;
 328     for (i = pairSearchList_.begin(); i != pairSearchList_.end(); ++i)
 329     {
 330         if (i->unique())
 331         {
 332             return *i;
 333         }
 334     }
 335     PairSearchImplPointer pairSearch(new AnalysisNeighborhoodPairSearchImpl(*this));
 336     pairSearchList_.push_back(pairSearch);
 337     return pairSearch;
 338 }
 339
 340 void AnalysisNeighborhoodSearchImpl::initGridCellNeighborList()
 341 {
 342     int   maxx, maxy, maxz;
 343     real  rvnorm;
 344
 345     /* Find the extent of the sphere in triclinic coordinates */
 346     maxz   = static_cast<int>(cutoff_ * recipcell_[ZZ][ZZ]) + 1;
 347     rvnorm = sqrt(sqr(recipcell_[YY][YY]) + sqr(recipcell_[ZZ][YY]));
 348     maxy   = static_cast<int>(cutoff_ * rvnorm) + 1;
 349     rvnorm = sqrt(sqr(recipcell_[XX][XX]) + sqr(recipcell_[YY][XX])
 350                   + sqr(recipcell_[ZZ][XX]));
 351     maxx   = static_cast<int>(cutoff_ * rvnorm) + 1;
 352
 353     /* Calculate the number of cells and reallocate if necessary */
 354     ngridnb_ = (2 * maxx + 1) * (2 * maxy + 1) * (2 * maxz + 1);
 355     if (gnboffs_nalloc_ < ngridnb_)
 356     {
 357         gnboffs_nalloc_ = ngridnb_;
 358         srenew(gnboffs_, gnboffs_nalloc_);
 359     }
 360
 361     /* Store the whole cube */
 362     /* TODO: Prune off corners that are not needed */
 363     int i = 0;
 364     for (int x = -maxx; x <= maxx; ++x)
 365     {
 366         for (int y = -maxy; y <= maxy; ++y)
 367         {
 368             for (int z = -maxz; z <= maxz; ++z)
 369             {
 370                 gnboffs_[i][XX] = x;
 371                 gnboffs_[i][YY] = y;
 372                 gnboffs_[i][ZZ] = z;
 373                 ++i;
 374             }
 375         }
 376     }
 377 }
 378
 379 bool AnalysisNeighborhoodSearchImpl::initGridCells(const t_pbc &pbc)
 380 {
 381     const real targetsize =
 382         pow(pbc.box[XX][XX] * pbc.box[YY][YY] * pbc.box[ZZ][ZZ]
 383             * 10 / nref_, static_cast<real>(1./3.));
 384
 385     int cellCount = 1;
 386     for (int dd = 0; dd < DIM; ++dd)
 387     {
 388         ncelldim_[dd] = static_cast<int>(pbc.box[dd][dd] / targetsize);
 389         cellCount    *= ncelldim_[dd];
 390         if (ncelldim_[dd] < 3)
 391         {
 392             return false;
 393         }
 394     }
 395     // Never decrease the size of the cell vector to avoid reallocating
 396     // memory for the nested vectors.  The actual size of the vector is not
 397     // used outside this function.
 398     if (cells_.size() < static_cast<size_t>(cellCount))
 399     {
 400         cells_.resize(cellCount);
 401     }
 402     for (int ci = 0; ci < cellCount; ++ci)
 403     {
 404         cells_[ci].clear();
 405     }
 406     return true;
 407 }
 408
 409 bool AnalysisNeighborhoodSearchImpl::initGrid(const t_pbc &pbc)
 410 {
 411     /* TODO: This check could be improved. */
 412     if (0.5*pbc.max_cutoff2 < cutoff2_)
 413     {
 414         return false;
 415     }
 416
 417     if (!initGridCells(pbc))
 418     {
 419         return false;
 420     }
 421
 422     bTric_ = TRICLINIC(pbc.box);
 423     if (bTric_)
 424     {
 425         for (int dd = 0; dd < DIM; ++dd)
 426         {
 427             svmul(1.0 / ncelldim_[dd], pbc.box[dd], cellbox_[dd]);
 428         }
 429         m_inv_ur0(cellbox_, recipcell_);
 430     }
 431     else
 432     {
 433         for (int dd = 0; dd < DIM; ++dd)
 434         {
 435             cellbox_[dd][dd]   = pbc.box[dd][dd] / ncelldim_[dd];
 436             recipcell_[dd][dd] = 1.0 / cellbox_[dd][dd];
 437         }
 438     }
 439     initGridCellNeighborList();
 440     return true;
 441 }
 442
 443 void AnalysisNeighborhoodSearchImpl::mapPointToGridCell(const rvec x,
 444                                                         ivec       cell,
 445                                                         rvec       xout) const
 446 {
 447     rvec xtmp;
 448     copy_rvec(x, xtmp);
 449     if (bTric_)
 450     {
 451         rvec tx;
 452         tmvmul_ur0(recipcell_, xtmp, tx);
 453         for (int dd = 0; dd < DIM; ++dd)
 454         {
 455             const int cellCount = ncelldim_[dd];
 456             int       cellIndex = static_cast<int>(floor(tx[dd]));
 457             while (cellIndex < 0)
 458             {
 459                 cellIndex += cellCount;
 460                 rvec_add(xtmp, pbc_.box[dd], xtmp);
 461             }
 462             while (cellIndex >= cellCount)
 463             {
 464                 cellIndex -= cellCount;
 465                 rvec_sub(xtmp, pbc_.box[dd], xtmp);
 466             }
 467             cell[dd] = cellIndex;
 468         }
 469     }
 470     else
 471     {
 472         for (int dd = 0; dd < DIM; ++dd)
 473         {
 474             const int cellCount = ncelldim_[dd];
 475             int       cellIndex = static_cast<int>(floor(xtmp[dd] * recipcell_[dd][dd]));
 476             while (cellIndex < 0)
 477             {
 478                 cellIndex += cellCount;
 479                 xtmp[dd]  += pbc_.box[dd][dd];
 480             }
 481             while (cellIndex >= cellCount)
 482             {
 483                 cellIndex -= cellCount;
 484                 xtmp[dd]  -= pbc_.box[dd][dd];
 485             }
 486             cell[dd] = cellIndex;
 487         }
 488     }
 489     copy_rvec(xtmp, xout);
 490 }
 491
 492 int AnalysisNeighborhoodSearchImpl::getGridCellIndex(const ivec cell) const
 493 {
 494     GMX_ASSERT(cell[XX] >= 0 && cell[XX] < ncelldim_[XX],
 495                "Grid cell X index out of range");
 496     GMX_ASSERT(cell[YY] >= 0 && cell[YY] < ncelldim_[YY],
 497                "Grid cell Y index out of range");
 498     GMX_ASSERT(cell[ZZ] >= 0 && cell[ZZ] < ncelldim_[ZZ],
 499                "Grid cell Z index out of range");
 500     return cell[XX] + cell[YY] * ncelldim_[XX]
 501            + cell[ZZ] * ncelldim_[XX] * ncelldim_[YY];
 502 }
 503
 504 void AnalysisNeighborhoodSearchImpl::addToGridCell(const ivec cell, int i)
 505 {
 506     const int ci = getGridCellIndex(cell);
 507     cells_[ci].push_back(i);
 508 }
 509
 510 void AnalysisNeighborhoodSearchImpl::init(
 511         AnalysisNeighborhood::SearchMode     mode,
 512         bool                                 bXY,
 513         const t_blocka                      *excls,
 514         const t_pbc                         *pbc,
 515         const AnalysisNeighborhoodPositions &positions)
 516 {
 517     GMX_RELEASE_ASSERT(positions.index_ == -1,
 518                        "Individual indexed positions not supported as reference");
 519     bXY_ = bXY;
 520     if (bXY_ && pbc->ePBC != epbcNONE)
 521     {
 522         if (pbc->ePBC != epbcXY && pbc->ePBC != epbcXYZ)
 523         {
 524             std::string message =
 525                 formatString("Computations in the XY plane are not supported with PBC type '%s'",
 526                              EPBC(pbc->ePBC));
 527             GMX_THROW(NotImplementedError(message));
 528         }
 529         if (std::fabs(pbc->box[ZZ][XX]) > GMX_REAL_EPS*pbc->box[ZZ][ZZ] ||
 530             std::fabs(pbc->box[ZZ][YY]) > GMX_REAL_EPS*pbc->box[ZZ][ZZ])
 531         {
 532             GMX_THROW(NotImplementedError("Computations in the XY plane are not supported when the last box vector is not parallel to the Z axis"));
 533         }
 534         set_pbc(&pbc_, epbcXY, const_cast<rvec *>(pbc->box));
 535     }
 536     else if (pbc != NULL)
 537     {
 538         pbc_  = *pbc;
 539     }
 540     else
 541     {
 542         pbc_.ePBC = epbcNONE;
 543     }
 544     nref_ = positions.count_;
 545     // TODO: Consider whether it would be possible to support grid searching in
 546     // more cases.
 547     if (mode == AnalysisNeighborhood::eSearchMode_Simple
 548         || pbc_.ePBC != epbcXYZ)
 549     {
 550         bGrid_ = false;
 551     }
 552     else if (bTryGrid_)
 553     {
 554         // TODO: Actually implement forcing eSearchMode_Grid
 555         bGrid_ = initGrid(pbc_);
 556     }
 557     if (bGrid_)
 558     {
 559         if (xref_nalloc_ < nref_)
 560         {
 561             srenew(xref_alloc_, nref_);
 562             xref_nalloc_ = nref_;
 563         }
 564         xref_ = xref_alloc_;
 565
 566         for (int i = 0; i < nref_; ++i)
 567         {
 568             ivec refcell;
 569             mapPointToGridCell(positions.x_[i], refcell, xref_alloc_[i]);
 570             addToGridCell(refcell, i);
 571         }
 572     }
 573     else
 574     {
 575         xref_ = positions.x_;
 576     }
 577     excls_           = excls;
 578     refExclusionIds_ = NULL;
 579     if (excls != NULL)
 580     {
 581         // TODO: Check that the IDs are ascending, or remove the limitation.
 582         refExclusionIds_ = positions.exclusionIds_;
 583         GMX_RELEASE_ASSERT(refExclusionIds_ != NULL,
 584                            "Exclusion IDs must be set for reference positions "
 585                            "when exclusions are enabled");
 586     }
 587 }
 588
 589 /********************************************************************
 590  * AnalysisNeighborhoodPairSearchImpl
 591  */
 592
 593 void AnalysisNeighborhoodPairSearchImpl::reset(int testIndex)
 594 {
 595     testIndex_ = testIndex;
 596     if (testIndex_ >= 0 && testIndex_ < static_cast<int>(testPositions_.size()))
 597     {
 598         copy_rvec(testPositions_[testIndex_], xtest_);
 599         if (search_.bGrid_)
 600         {
 601             search_.mapPointToGridCell(testPositions_[testIndex], testcell_, xtest_);
 602         }
 603         else
 604         {
 605             copy_rvec(testPositions_[testIndex_], xtest_);
 606         }
 607         if (search_.excls_ != NULL)
 608         {
 609             const int exclIndex  = testExclusionIds_[testIndex];
 610             if (exclIndex < search_.excls_->nr)
 611             {
 612                 const int startIndex = search_.excls_->index[exclIndex];
 613                 nexcl_ = search_.excls_->index[exclIndex + 1] - startIndex;
 614                 excl_  = &search_.excls_->a[startIndex];
 615             }
 616             else
 617             {
 618                 nexcl_ = 0;
 619                 excl_  = NULL;
 620             }
 621         }
 622     }
 623     previ_     = -1;
 624     prevr2_    = 0.0;
 625     exclind_   = 0;
 626     prevnbi_   = 0;
 627     prevcai_   = -1;
 628 }
 629
 630 void AnalysisNeighborhoodPairSearchImpl::nextTestPosition()
 631 {
 632     if (testIndex_ < static_cast<int>(testPositions_.size()))
 633     {
 634         ++testIndex_;
 635         reset(testIndex_);
 636     }
 637 }
 638
 639 bool AnalysisNeighborhoodPairSearchImpl::isExcluded(int j)
 640 {
 641     if (exclind_ < nexcl_)
 642     {
 643         const int refId = search_.refExclusionIds_[j];
 644         while (exclind_ < nexcl_ && excl_[exclind_] < refId)
 645         {
 646             ++exclind_;
 647         }
 648         if (exclind_ < nexcl_ && refId == excl_[exclind_])
 649         {
 650             ++exclind_;
 651             return true;
 652         }
 653     }
 654     return false;
 655 }
 656
 657 void AnalysisNeighborhoodPairSearchImpl::startSearch(
 658         const AnalysisNeighborhoodPositions &positions)
 659 {
 660     testExclusionIds_ = positions.exclusionIds_;
 661     GMX_RELEASE_ASSERT(search_.excls_ == NULL || testExclusionIds_ != NULL,
 662                        "Exclusion IDs must be set when exclusions are enabled");
 663     if (positions.index_ < 0)
 664     {
 665         testPositions_ = constArrayRefFromArray<rvec>(positions.x_, positions.count_);
 666         reset(0);
 667     }
 668     else
 669     {
 670         // Somewhat of a hack: setup the array such that only the last position
 671         // will be used.
 672         testPositions_ = constArrayRefFromArray<rvec>(positions.x_, positions.index_ + 1);
 673         reset(positions.index_);
 674     }
 675 }
 676
 677 template <class Action>
 678 bool AnalysisNeighborhoodPairSearchImpl::searchNext(Action action)
 679 {
 680     while (testIndex_ < static_cast<int>(testPositions_.size()))
 681     {
 682         if (search_.bGrid_)
 683         {
 684             GMX_RELEASE_ASSERT(!search_.bXY_, "Grid-based XY searches not implemented");
 685
 686             int nbi = prevnbi_;
 687             int cai = prevcai_ + 1;
 688
 689             for (; nbi < search_.ngridnb_; ++nbi)
 690             {
 691                 ivec cell;
 692
 693                 ivec_add(testcell_, search_.gnboffs_[nbi], cell);
 694                 cell[XX] = (cell[XX] + search_.ncelldim_[XX]) % search_.ncelldim_[XX];
 695                 cell[YY] = (cell[YY] + search_.ncelldim_[YY]) % search_.ncelldim_[YY];
 696                 cell[ZZ] = (cell[ZZ] + search_.ncelldim_[ZZ]) % search_.ncelldim_[ZZ];
 697
 698                 const int ci       = search_.getGridCellIndex(cell);
 699                 const int cellSize = static_cast<int>(search_.cells_[ci].size());
 700                 /* TODO: Calculate the required PBC shift outside the inner loop */
 701                 for (; cai < cellSize; ++cai)
 702                 {
 703                     const int i = search_.cells_[ci][cai];
 704                     if (isExcluded(i))
 705                     {
 706                         continue;
 707                     }
 708                     rvec       dx;
 709                     pbc_dx_aiuc(&search_.pbc_, xtest_, search_.xref_[i], dx);
 710                     const real r2 = norm2(dx);
 711                     if (r2 <= search_.cutoff2_)
 712                     {
 713                         if (action(i, r2))
 714                         {
 715                             prevnbi_ = nbi;
 716                             prevcai_ = cai;
 717                             previ_   = i;
 718                             prevr2_  = r2;
 719                             return true;
 720                         }
 721                     }
 722                 }
 723                 exclind_ = 0;
 724                 cai      = 0;
 725             }
 726         }
 727         else
 728         {
 729             for (int i = previ_ + 1; i < search_.nref_; ++i)
 730             {
 731                 if (isExcluded(i))
 732                 {
 733                     continue;
 734                 }
 735                 rvec dx;
 736                 if (search_.pbc_.ePBC != epbcNONE)
 737                 {
 738                     pbc_dx(&search_.pbc_, xtest_, search_.xref_[i], dx);
 739                 }
 740                 else
 741                 {
 742                     rvec_sub(xtest_, search_.xref_[i], dx);
 743                 }
 744                 const real r2
 745                     = search_.bXY_
 746                         ? dx[XX]*dx[XX] + dx[YY]*dx[YY]
 747                         : norm2(dx);
 748                 if (r2 <= search_.cutoff2_)
 749                 {
 750                     if (action(i, r2))
 751                     {
 752                         previ_  = i;
 753                         prevr2_ = r2;
 754                         return true;
 755                     }
 756                 }
 757             }
 758         }
 759         nextTestPosition();
 760     }
 761     return false;
 762 }
 763
 764 void AnalysisNeighborhoodPairSearchImpl::initFoundPair(
 765         AnalysisNeighborhoodPair *pair) const
 766 {
 767     if (previ_ < 0)
 768     {
 769         *pair = AnalysisNeighborhoodPair();
 770     }
 771     else
 772     {
 773         *pair = AnalysisNeighborhoodPair(previ_, testIndex_, prevr2_);
 774     }
 775 }
 776
 777 }   // namespace internal
 778
 779 namespace
 780 {
 781
 782 /*! \brief
 783  * Search action to find the next neighbor.
 784  *
 785  * Used as the action for AnalysisNeighborhoodPairSearchImpl::searchNext() to
 786  * find the next neighbor.
 787  *
 788  * Simply breaks the loop on the first found neighbor.
 789  */
 790 bool withinAction(int /*i*/, real /*r2*/)
 791 {
 792     return true;
 793 }
 794
 795 /*! \brief
 796  * Search action find the minimum distance.
 797  *
 798  * Used as the action for AnalysisNeighborhoodPairSearchImpl::searchNext() to
 799  * find the nearest neighbor.
 800  *
 801  * With this action, AnalysisNeighborhoodPairSearchImpl::searchNext() always
 802  * returns false, and the output is put into the variables passed by pointer
 803  * into the constructor.  If no neighbors are found, the output variables are
 804  * not modified, i.e., the caller must initialize them.
 805  */
 806 class MindistAction
 807 {
 808     public:
 809         /*! \brief
 810          * Initializes the action with given output locations.
 811          *
 812          * \param[out] closestPoint Index of the closest reference location.
 813          * \param[out] minDist2     Minimum distance squared.
 814          *
 815          * The constructor call does not modify the pointed values, but only
 816          * stores the pointers for later use.
 817          * See the class description for additional semantics.
 818          */
 819         MindistAction(int *closestPoint, real *minDist2)
 820             : closestPoint_(*closestPoint), minDist2_(*minDist2)
 821         {
 822         }
 823
 824         //! Processes a neighbor to find the nearest point.
 825         bool operator()(int i, real r2)
 826         {
 827             if (r2 < minDist2_)
 828             {
 829                 closestPoint_ = i;
 830                 minDist2_     = r2;
 831             }
 832             return false;
 833         }
 834
 835     private:
 836         int     &closestPoint_;
 837         real    &minDist2_;
 838
 839         GMX_DISALLOW_ASSIGN(MindistAction);
 840 };
 841
 842 }   // namespace
 843
 844 /********************************************************************
 845  * AnalysisNeighborhood::Impl
 846  */
 847
 848 class AnalysisNeighborhood::Impl
 849 {
 850     public:
 851         typedef AnalysisNeighborhoodSearch::ImplPointer SearchImplPointer;
 852         typedef std::vector<SearchImplPointer> SearchList;
 853
 854         Impl() : cutoff_(0), excls_(NULL), mode_(eSearchMode_Automatic), bXY_(false)
 855         {
 856         }
 857         ~Impl()
 858         {
 859             SearchList::const_iterator i;
 860             for (i = searchList_.begin(); i != searchList_.end(); ++i)
 861             {
 862                 GMX_RELEASE_ASSERT(i->unique(),
 863                                    "Dangling AnalysisNeighborhoodSearch reference");
 864             }
 865         }
 866
 867         SearchImplPointer getSearch();
 868
 869         tMPI::mutex             createSearchMutex_;
 870         SearchList              searchList_;
 871         real                    cutoff_;
 872         const t_blocka         *excls_;
 873         SearchMode              mode_;
 874         bool                    bXY_;
 875 };
 876
 877 AnalysisNeighborhood::Impl::SearchImplPointer
 878 AnalysisNeighborhood::Impl::getSearch()
 879 {
 880     tMPI::lock_guard<tMPI::mutex> lock(createSearchMutex_);
 881     // TODO: Consider whether this needs to/can be faster, e.g., by keeping a
 882     // separate pool of unused search objects.
 883     SearchList::const_iterator i;
 884     for (i = searchList_.begin(); i != searchList_.end(); ++i)
 885     {
 886         if (i->unique())
 887         {
 888             return *i;
 889         }
 890     }
 891     SearchImplPointer search(new internal::AnalysisNeighborhoodSearchImpl(cutoff_));
 892     searchList_.push_back(search);
 893     return search;
 894 }
 895
 896 /********************************************************************
 897  * AnalysisNeighborhood
 898  */
 899
 900 AnalysisNeighborhood::AnalysisNeighborhood()
 901     : impl_(new Impl)
 902 {
 903 }
 904
 905 AnalysisNeighborhood::~AnalysisNeighborhood()
 906 {
 907 }
 908
 909 void AnalysisNeighborhood::setCutoff(real cutoff)
 910 {
 911     GMX_RELEASE_ASSERT(impl_->searchList_.empty(),
 912                        "Changing the cutoff after initSearch() not currently supported");
 913     impl_->cutoff_ = cutoff;
 914 }
 915
 916 void AnalysisNeighborhood::setXYMode(bool bXY)
 917 {
 918     impl_->bXY_ = bXY;
 919 }
 920
 921 void AnalysisNeighborhood::setTopologyExclusions(const t_blocka *excls)
 922 {
 923     GMX_RELEASE_ASSERT(impl_->searchList_.empty(),
 924                        "Changing the exclusions after initSearch() not currently supported");
 925     impl_->excls_ = excls;
 926 }
 927
 928 void AnalysisNeighborhood::setMode(SearchMode mode)
 929 {
 930     impl_->mode_ = mode;
 931 }
 932
 933 AnalysisNeighborhood::SearchMode AnalysisNeighborhood::mode() const
 934 {
 935     return impl_->mode_;
 936 }
 937
 938 AnalysisNeighborhoodSearch
 939 AnalysisNeighborhood::initSearch(const t_pbc                         *pbc,
 940                                  const AnalysisNeighborhoodPositions &positions)
 941 {
 942     Impl::SearchImplPointer search(impl_->getSearch());
 943     search->init(mode(), impl_->bXY_, impl_->excls_, pbc, positions);
 944     return AnalysisNeighborhoodSearch(search);
 945 }
 946
 947 /********************************************************************
 948  * AnalysisNeighborhoodSearch
 949  */
 950
 951 AnalysisNeighborhoodSearch::AnalysisNeighborhoodSearch()
 952 {
 953 }
 954
 955 AnalysisNeighborhoodSearch::AnalysisNeighborhoodSearch(const ImplPointer &impl)
 956     : impl_(impl)
 957 {
 958 }
 959
 960 void AnalysisNeighborhoodSearch::reset()
 961 {
 962     impl_.reset();
 963 }
 964
 965 AnalysisNeighborhood::SearchMode AnalysisNeighborhoodSearch::mode() const
 966 {
 967     GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
 968     return (impl_->usesGridSearch()
 969             ? AnalysisNeighborhood::eSearchMode_Grid
 970             : AnalysisNeighborhood::eSearchMode_Simple);
 971 }
 972
 973 bool AnalysisNeighborhoodSearch::isWithin(
 974         const AnalysisNeighborhoodPositions &positions) const
 975 {
 976     GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
 977     internal::AnalysisNeighborhoodPairSearchImpl pairSearch(*impl_);
 978     pairSearch.startSearch(positions);
 979     return pairSearch.searchNext(&withinAction);
 980 }
 981
 982 real AnalysisNeighborhoodSearch::minimumDistance(
 983         const AnalysisNeighborhoodPositions &positions) const
 984 {
 985     GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
 986     internal::AnalysisNeighborhoodPairSearchImpl pairSearch(*impl_);
 987     pairSearch.startSearch(positions);
 988     real          minDist2     = impl_->cutoffSquared();
 989     int           closestPoint = -1;
 990     MindistAction action(&closestPoint, &minDist2);
 991     (void)pairSearch.searchNext(action);
 992     return sqrt(minDist2);
 993 }
 994
 995 AnalysisNeighborhoodPair
 996 AnalysisNeighborhoodSearch::nearestPoint(
 997         const AnalysisNeighborhoodPositions &positions) const
 998 {
 999     GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
1000     internal::AnalysisNeighborhoodPairSearchImpl pairSearch(*impl_);
1001     pairSearch.startSearch(positions);
1002     real          minDist2     = impl_->cutoffSquared();
1003     int           closestPoint = -1;
1004     MindistAction action(&closestPoint, &minDist2);
1005     (void)pairSearch.searchNext(action);
1006     return AnalysisNeighborhoodPair(closestPoint, 0, minDist2);
1007 }
1008
1009 AnalysisNeighborhoodPairSearch
1010 AnalysisNeighborhoodSearch::startPairSearch(
1011         const AnalysisNeighborhoodPositions &positions) const
1012 {
1013     GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
1014     Impl::PairSearchImplPointer pairSearch(impl_->getPairSearch());
1015     pairSearch->startSearch(positions);
1016     return AnalysisNeighborhoodPairSearch(pairSearch);
1017 }
1018
1019 /********************************************************************
1020  * AnalysisNeighborhoodPairSearch
1021  */
1022
1023 AnalysisNeighborhoodPairSearch::AnalysisNeighborhoodPairSearch(
1024         const ImplPointer &impl)
1025     : impl_(impl)
1026 {
1027 }
1028
1029 bool AnalysisNeighborhoodPairSearch::findNextPair(AnalysisNeighborhoodPair *pair)
1030 {
1031     bool bFound = impl_->searchNext(&withinAction);
1032     impl_->initFoundPair(pair);
1033     return bFound;
1034 }
1035
1036 void AnalysisNeighborhoodPairSearch::skipRemainingPairsForTestPosition()
1037 {
1038     impl_->nextTestPosition();
1039 }
1040
1041 } // namespace gmx