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,2015,2016, 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  * High-level overview of the algorithm is at \ref page_analysisnbsearch.
  40  *
  41  * \todo
  42  * The grid implementation could still be optimized in several different ways:
  43  *   - Pruning grid cells from the search list if they are completely outside
  44  *     the sphere that is being considered.
  45  *   - A better heuristic could be added for falling back to simple loops for a
  46  *     small number of reference particles.
  47  *   - A better heuristic for selecting the grid size.
  48  *   - A multi-level grid implementation could be used to be able to use small
  49  *     grids for short cutoffs with very inhomogeneous particle distributions
  50  *     without a memory cost.
  51  *
  52  * \author Teemu Murtola <teemu.murtola@gmail.com>
  53  * \ingroup module_selection
  54  */
  55 #include "gmxpre.h"
  56
  57 #include "nbsearch.h"
  58
  59 #include <cmath>
  60 #include <cstring>
  61
  62 #include <algorithm>
  63 #include <vector>
  64
  65 #include "gromacs/math/functions.h"
  66 #include "gromacs/math/vec.h"
  67 #include "gromacs/pbcutil/pbc.h"
  68 #include "gromacs/selection/position.h"
  69 #include "gromacs/topology/block.h"
  70 #include "gromacs/utility/arrayref.h"
  71 #include "gromacs/utility/exceptions.h"
  72 #include "gromacs/utility/gmxassert.h"
  73 #include "gromacs/utility/mutex.h"
  74 #include "gromacs/utility/stringutil.h"
  75
  76 namespace gmx
  77 {
  78
  79 namespace
  80 {
  81
  82 /*! \brief
  83  * Computes the bounding box for a set of positions.
  84  *
  85  * \param[in]  posCount Number of positions in \p x.
  86  * \param[in]  x        Positions to compute the bounding box for.
  87  * \param[out] origin   Origin of the bounding box.
  88  * \param[out] size     Size of the bounding box.
  89  */
  90 void computeBoundingBox(int posCount, const rvec x[], rvec origin, rvec size)
  91 {
  92     rvec maxBound;
  93     copy_rvec(x[0], origin);
  94     copy_rvec(x[0], maxBound);
  95     for (int i = 1; i < posCount; ++i)
  96     {
  97         for (int d = 0; d < DIM; ++d)
  98         {
  99             if (origin[d] > x[i][d])
 100             {
 101                 origin[d] = x[i][d];
 102             }
 103             if (maxBound[d] < x[i][d])
 104             {
 105                 maxBound[d] = x[i][d];
 106             }
 107         }
 108     }
 109     rvec_sub(maxBound, origin, size);
 110 }
 111
 112 }   // namespace
 113
 114 namespace internal
 115 {
 116
 117 /********************************************************************
 118  * Implementation class declarations
 119  */
 120
 121 class AnalysisNeighborhoodSearchImpl
 122 {
 123     public:
 124         typedef AnalysisNeighborhoodPairSearch::ImplPointer
 125             PairSearchImplPointer;
 126         typedef std::vector<PairSearchImplPointer> PairSearchList;
 127         typedef std::vector<std::vector<int> > CellList;
 128
 129         explicit AnalysisNeighborhoodSearchImpl(real cutoff);
 130         ~AnalysisNeighborhoodSearchImpl();
 131
 132         /*! \brief
 133          * Initializes the search with a given box and reference positions.
 134          *
 135          * \param[in] mode            Search mode to use.
 136          * \param[in] bXY             Whether to use 2D searching.
 137          * \param[in] excls           Exclusions.
 138          * \param[in] pbc             PBC information.
 139          * \param[in] positions       Set of reference positions.
 140          */
 141         void init(AnalysisNeighborhood::SearchMode     mode,
 142                   bool                                 bXY,
 143                   const t_blocka                      *excls,
 144                   const t_pbc                         *pbc,
 145                   const AnalysisNeighborhoodPositions &positions);
 146         PairSearchImplPointer getPairSearch();
 147
 148         real cutoffSquared() const { return cutoff2_; }
 149         bool usesGridSearch() const { return bGrid_; }
 150
 151     private:
 152         /*! \brief
 153          * Checks the efficiency and possibility of doing grid-based searching.
 154          *
 155          * \param[in] bForce  If `true`, grid search will be forced if possible.
 156          * \returns   `false` if grid search is not suitable.
 157          */
 158         bool checkGridSearchEfficiency(bool bForce);
 159         /*! \brief
 160          * Determines a suitable grid size and sets up the cells.
 161          *
 162          * \param[in] box          Box vectors (should not have zero vectors).
 163          * \param[in] bSingleCell  If `true`, the corresponding dimension will
 164          *     be forced to use a single cell.
 165          * \param[in] posCount     Number of positions that will be put on the
 166          *     grid.
 167          * \returns   `false` if grid search is not suitable.
 168          */
 169         bool initGridCells(const matrix box, bool bSingleCell[DIM],
 170                            int posCount);
 171         /*! \brief
 172          * Sets ua a search grid for a given box.
 173          *
 174          * \param[in] pbc      Information about the box.
 175          * \param[in] posCount Number of positions in \p x.
 176          * \param[in] x        Reference positions that will be put on the grid.
 177          * \param[in] bForce   If `true`, grid searching will be used if at all
 178          *     possible, even if a simple search might give better performance.
 179          * \returns   `false` if grid search is not suitable.
 180          */
 181         bool initGrid(const t_pbc &pbc, int posCount, const rvec x[], bool bForce);
 182         /*! \brief
 183          * Maps a point into a grid cell.
 184          *
 185          * \param[in]  x    Point to map.
 186          * \param[out] cell Fractional cell coordinates of \p x on the grid.
 187          * \param[out] xout Coordinates to use.
 188          *
 189          * \p xout will be within the rectangular unit cell in dimensions where
 190          * the grid is periodic.  For other dimensions, both \p xout and
 191          * \p cell can be outside the grid/unit cell.
 192          */
 193         void mapPointToGridCell(const rvec x, rvec cell, rvec xout) const;
 194         /*! \brief
 195          * Calculates linear index of a grid cell.
 196          *
 197          * \param[in]  cell Cell indices (must be within the grid).
 198          * \returns    Linear index of \p cell.
 199          */
 200         int getGridCellIndex(const ivec cell) const;
 201         /*! \brief
 202          * Adds an index into a grid cell.
 203          *
 204          * \param[in]  cell Fractional cell coordinates into which \p i should
 205          *     be added.
 206          * \param[in]  i    Index to add.
 207          *
 208          * \p cell should satisfy the conditions that \p mapPointToGridCell()
 209          * produces.
 210          */
 211         void addToGridCell(const rvec cell, int i);
 212         /*! \brief
 213          * Initializes a cell pair loop for a dimension.
 214          *
 215          * \param[in]     centerCell Fractional cell coordiates of the particle
 216          *     for which pairs are being searched.
 217          * \param[in,out] cell       Current/initial cell to loop over.
 218          * \param[in,out] upperBound Last cell to loop over.
 219          * \param[in]     dim        Dimension to initialize in this call.
 220          *
 221          * Initializes `cell[dim]` and `upperBound[dim]` for looping over
 222          * neighbors of a particle at position given by \p centerCell.
 223          * If 'dim != ZZ`, `cell[d]` (`d > dim`) set the plane/row of cells
 224          * for which the loop is initialized.  The loop should then go from
 225          * `cell[dim]` until `upperBound[dim]`, inclusive.
 226          * `cell[d]` with `d < dim` or `upperBound[d]` with `d != dim` are not
 227          * modified by this function.
 228          *
 229          * `cell` and `upperBound` may be outside the grid for periodic
 230          * dimensions and need to be shifted separately: to simplify the
 231          * looping, the range is always (roughly) symmetric around the value in
 232          * `centerCell`.
 233          */
 234         void initCellRange(const rvec centerCell, ivec cell,
 235                            ivec upperBound, int dim) const;
 236         /*! \brief
 237          * Advances cell pair loop to the next cell.
 238          *
 239          * \param[in]     centerCell Fractional cell coordiates of the particle
 240          *     for which pairs are being searched.
 241          * \param[in,out] cell       Current (in)/next (out) cell in the loop.
 242          * \param[in,out] upperBound Last cell in the loop for each dimension.
 243          */
 244         bool nextCell(const rvec centerCell, ivec cell, ivec upperBound) const;
 245         /*! \brief
 246          * Calculates the index and shift of a grid cell during looping.
 247          *
 248          * \param[in]  cell       Unshifted cell index.
 249          * \param[out] shift      Shift to apply to get the periodic distance
 250          *     for distances between the cells.
 251          * \returns    Grid cell index corresponding to `cell`.
 252          */
 253         int shiftCell(const ivec cell, rvec shift) const;
 254
 255         //! Whether to try grid searching.
 256         bool                    bTryGrid_;
 257         //! The cutoff.
 258         real                    cutoff_;
 259         //! The cutoff squared.
 260         real                    cutoff2_;
 261         //! Whether to do searching in XY plane only.
 262         bool                    bXY_;
 263
 264         //! Number of reference points for the current frame.
 265         int                     nref_;
 266         //! Reference point positions.
 267         const rvec             *xref_;
 268         //! Reference position exclusion IDs.
 269         const int              *refExclusionIds_;
 270         //! Reference position indices (NULL if no indices).
 271         const int              *refIndices_;
 272         //! Exclusions.
 273         const t_blocka         *excls_;
 274         //! PBC data.
 275         t_pbc                   pbc_;
 276
 277         //! Whether grid searching is actually used for the current positions.
 278         bool                    bGrid_;
 279         //! false if the box is rectangular.
 280         bool                    bTric_;
 281         //! Whether the grid is periodic in a dimension.
 282         bool                    bGridPBC_[DIM];
 283         //! Array for storing in-unit-cell reference positions.
 284         std::vector<RVec>       xrefAlloc_;
 285         //! Origin of the grid (zero for periodic dimensions).
 286         rvec                    gridOrigin_;
 287         //! Size of a single grid cell.
 288         rvec                    cellSize_;
 289         //! Inverse of \p cellSize_. Zero for dimensions where grid is not used.
 290         rvec                    invCellSize_;
 291         /*! \brief
 292          * Shift in cell coordinates (for triclinic boxes) in X when crossing
 293          * the Z periodic boundary.
 294          */
 295         real                    cellShiftZX_;
 296         /*! \brief
 297          * Shift in cell coordinates (for triclinic boxes) in Y when crossing
 298          * the Z periodic boundary.
 299          */
 300         real                    cellShiftZY_;
 301         /*! \brief
 302          * Shift in cell coordinates (for triclinic boxes) in X when crossing
 303          * the Y periodic boundary.
 304          */
 305         real                    cellShiftYX_;
 306         //! Number of cells along each dimension.
 307         ivec                    ncelldim_;
 308         //! Data structure to hold the grid cell contents.
 309         CellList                cells_;
 310
 311         Mutex                   createPairSearchMutex_;
 312         PairSearchList          pairSearchList_;
 313
 314         friend class AnalysisNeighborhoodPairSearchImpl;
 315
 316         GMX_DISALLOW_COPY_AND_ASSIGN(AnalysisNeighborhoodSearchImpl);
 317 };
 318
 319 class AnalysisNeighborhoodPairSearchImpl
 320 {
 321     public:
 322         explicit AnalysisNeighborhoodPairSearchImpl(const AnalysisNeighborhoodSearchImpl &search)
 323             : search_(search)
 324         {
 325             testPosCount_     = 0;
 326             testPositions_    = NULL;
 327             testExclusionIds_ = NULL;
 328             testIndices_      = NULL;
 329             nexcl_            = 0;
 330             excl_             = NULL;
 331             clear_rvec(xtest_);
 332             clear_rvec(testcell_);
 333             clear_ivec(currCell_);
 334             clear_ivec(cellBound_);
 335             reset(-1);
 336         }
 337
 338         //! Initializes a search to find reference positions neighboring \p x.
 339         void startSearch(const AnalysisNeighborhoodPositions &positions);
 340         //! Searches for the next neighbor.
 341         template <class Action>
 342         bool searchNext(Action action);
 343         //! Initializes a pair representing the pair found by searchNext().
 344         void initFoundPair(AnalysisNeighborhoodPair *pair) const;
 345         //! Advances to the next test position, skipping any remaining pairs.
 346         void nextTestPosition();
 347
 348     private:
 349         //! Clears the loop indices.
 350         void reset(int testIndex);
 351         //! Checks whether a reference positiong should be excluded.
 352         bool isExcluded(int j);
 353
 354         //! Parent search object.
 355         const AnalysisNeighborhoodSearchImpl   &search_;
 356         //! Number of test positions.
 357         int                                     testPosCount_;
 358         //! Reference to the test positions.
 359         const rvec                             *testPositions_;
 360         //! Reference to the test exclusion indices.
 361         const int                              *testExclusionIds_;
 362         //! Reference to the test position indices.
 363         const int                              *testIndices_;
 364         //! Number of excluded reference positions for current test particle.
 365         int                                     nexcl_;
 366         //! Exclusions for current test particle.
 367         const int                              *excl_;
 368         //! Index of the currently active test position in \p testPositions_.
 369         int                                     testIndex_;
 370         //! Stores test position during a pair loop.
 371         rvec                                    xtest_;
 372         //! Stores the previous returned position during a pair loop.
 373         int                                     previ_;
 374         //! Stores the pair distance corresponding to previ_;
 375         real                                    prevr2_;
 376         //! Stores the shortest distance vector corresponding to previ_;
 377         rvec                                    prevdx_;
 378         //! Stores the current exclusion index during loops.
 379         int                                     exclind_;
 380         //! Stores the fractional test particle cell location during loops.
 381         rvec                                    testcell_;
 382         //! Stores the current cell during pair loops.
 383         ivec                                    currCell_;
 384         //! Stores the current loop upper bounds for each dimension during pair loops.
 385         ivec                                    cellBound_;
 386         //! Stores the index within the current cell during pair loops.
 387         int                                     prevcai_;
 388
 389         GMX_DISALLOW_COPY_AND_ASSIGN(AnalysisNeighborhoodPairSearchImpl);
 390 };
 391
 392 /********************************************************************
 393  * AnalysisNeighborhoodSearchImpl
 394  */
 395
 396 AnalysisNeighborhoodSearchImpl::AnalysisNeighborhoodSearchImpl(real cutoff)
 397 {
 398     bTryGrid_       = true;
 399     cutoff_         = cutoff;
 400     if (cutoff_ <= 0)
 401     {
 402         cutoff_     = cutoff2_ = GMX_REAL_MAX;
 403         bTryGrid_   = false;
 404     }
 405     else
 406     {
 407         cutoff2_        = gmx::square(cutoff_);
 408     }
 409     bXY_             = false;
 410     nref_            = 0;
 411     xref_            = NULL;
 412     refExclusionIds_ = NULL;
 413     refIndices_      = NULL;
 414     std::memset(&pbc_, 0, sizeof(pbc_));
 415
 416     bGrid_          = false;
 417     bTric_          = false;
 418     bGridPBC_[XX]   = true;
 419     bGridPBC_[YY]   = true;
 420     bGridPBC_[ZZ]   = true;
 421
 422     clear_rvec(gridOrigin_);
 423     clear_rvec(cellSize_);
 424     clear_rvec(invCellSize_);
 425     clear_ivec(ncelldim_);
 426 }
 427
 428 AnalysisNeighborhoodSearchImpl::~AnalysisNeighborhoodSearchImpl()
 429 {
 430     PairSearchList::const_iterator i;
 431     for (i = pairSearchList_.begin(); i != pairSearchList_.end(); ++i)
 432     {
 433         GMX_RELEASE_ASSERT(i->unique(),
 434                            "Dangling AnalysisNeighborhoodPairSearch reference");
 435     }
 436 }
 437
 438 AnalysisNeighborhoodSearchImpl::PairSearchImplPointer
 439 AnalysisNeighborhoodSearchImpl::getPairSearch()
 440 {
 441     lock_guard<Mutex> lock(createPairSearchMutex_);
 442     // TODO: Consider whether this needs to/can be faster, e.g., by keeping a
 443     // separate pool of unused search objects.
 444     PairSearchList::const_iterator i;
 445     for (i = pairSearchList_.begin(); i != pairSearchList_.end(); ++i)
 446     {
 447         if (i->unique())
 448         {
 449             return *i;
 450         }
 451     }
 452     PairSearchImplPointer pairSearch(new AnalysisNeighborhoodPairSearchImpl(*this));
 453     pairSearchList_.push_back(pairSearch);
 454     return pairSearch;
 455 }
 456
 457 bool AnalysisNeighborhoodSearchImpl::checkGridSearchEfficiency(bool bForce)
 458 {
 459     // Find the extent of the sphere in cells.
 460     ivec  range;
 461     for (int dd = 0; dd < DIM; ++dd)
 462     {
 463         range[dd] = static_cast<int>(ceil(cutoff_ * invCellSize_[dd]));
 464     }
 465
 466     // Calculate the fraction of cell pairs that need to be searched,
 467     // and check that the cutoff is not too large for periodic dimensions.
 468     real coveredCells = 1.0;
 469     for (int dd = 0; dd < DIM; ++dd)
 470     {
 471         const int cellCount    = ncelldim_[dd];
 472         const int coveredCount = 2 * range[dd] + 1;
 473         if (bGridPBC_[dd])
 474         {
 475             if (coveredCount > cellCount)
 476             {
 477                 // Cutoff is too close to half the box size for grid searching
 478                 // (it is not possible to find a single shift for every pair of
 479                 // grid cells).
 480                 return false;
 481             }
 482             coveredCells *= coveredCount;
 483         }
 484         else
 485         {
 486             if (range[dd] >= cellCount - 1)
 487             {
 488                 range[dd]     = cellCount - 1;
 489                 coveredCells *= cellCount;
 490             }
 491             else if (coveredCount > cellCount)
 492             {
 493                 // The sum of range+1, range+2, ..., range+N/2, ... range+1.
 494                 coveredCells *= range[dd] +
 495                     static_cast<real>((cellCount + 1)/2 * (cellCount/2 + 1)) / cellCount;
 496             }
 497             else
 498             {
 499                 // The sum of range+1, ..., 2*range+1, ..., 2*range+1, ... range+1.
 500                 coveredCells *= coveredCount -
 501                     static_cast<real>(range[dd] * (range[dd] + 1)) / cellCount;
 502             }
 503         }
 504     }
 505     // Magic constant that would need tuning for optimal performance:
 506     // Don't do grid searching if nearly all cell pairs would anyways need to
 507     // be looped through.
 508     const int totalCellCount = ncelldim_[XX] * ncelldim_[YY] * ncelldim_[ZZ];
 509     if (!bForce && coveredCells >= 0.5 * totalCellCount)
 510     {
 511         return false;
 512     }
 513     return true;
 514 }
 515
 516 bool AnalysisNeighborhoodSearchImpl::initGridCells(
 517         const matrix box, bool bSingleCell[DIM], int posCount)
 518 {
 519     // Determine the size of cubes where there are on average 10 positions.
 520     // The loop takes care of cases where some of the box edges are shorter
 521     // than the the desired cube size; in such cases, a single grid cell is
 522     // used in these dimensions, and the cube size is determined only from the
 523     // larger box vectors.  Such boxes should be rare, but the bounding box
 524     // approach can result in very flat boxes with certain types of selections
 525     // (e.g., for interfacial systems or for small number of atoms).
 526     real targetsize   = 0.0;
 527     int  prevDimCount = 4;
 528     while (true)
 529     {
 530         real volume   = 1.0;
 531         int  dimCount = 3;
 532         for (int dd = 0; dd < DIM; ++dd)
 533         {
 534             const real boxSize = box[dd][dd];
 535             if (boxSize < targetsize)
 536             {
 537                 bSingleCell[dd] = true;
 538                 if (bGridPBC_[dd])
 539                 {
 540                     return false;
 541                 }
 542             }
 543             if (bSingleCell[dd])
 544             {
 545                 --dimCount;
 546             }
 547             else
 548             {
 549                 volume *= boxSize;
 550             }
 551         }
 552         if (dimCount == 0 || dimCount == prevDimCount)
 553         {
 554             break;
 555         }
 556         targetsize   = pow(volume * 10 / posCount, static_cast<real>(1./dimCount));
 557         prevDimCount = dimCount;
 558     }
 559
 560     int totalCellCount = 1;
 561     for (int dd = 0; dd < DIM; ++dd)
 562     {
 563         int cellCount;
 564         if (bSingleCell[dd])
 565         {
 566             cellCount = 1;
 567         }
 568         else
 569         {
 570             cellCount = std::max(1, static_cast<int>(box[dd][dd] / targetsize));
 571             // TODO: If the cell count is one or two, it would be better to
 572             // just fall back to bSingleCell[dd] = true, and leave the rest to
 573             // the efficiency check later.
 574             if (bGridPBC_[dd] && cellCount < 3)
 575             {
 576                 return false;
 577             }
 578         }
 579         totalCellCount *= cellCount;
 580         ncelldim_[dd]   = cellCount;
 581     }
 582     if (totalCellCount <= 3)
 583     {
 584         return false;
 585     }
 586     // Never decrease the size of the cell vector to avoid reallocating
 587     // memory for the nested vectors.  The actual size of the vector is not
 588     // used outside this function.
 589     if (cells_.size() < static_cast<size_t>(totalCellCount))
 590     {
 591         cells_.resize(totalCellCount);
 592     }
 593     for (int ci = 0; ci < totalCellCount; ++ci)
 594     {
 595         cells_[ci].clear();
 596     }
 597     return true;
 598 }
 599
 600 bool AnalysisNeighborhoodSearchImpl::initGrid(
 601         const t_pbc &pbc, int posCount, const rvec x[], bool bForce)
 602 {
 603     if (posCount == 0)
 604     {
 605         return false;
 606     }
 607
 608     switch (pbc.ePBC)
 609     {
 610         case epbcNONE:
 611             bGridPBC_[XX] = false;
 612             bGridPBC_[YY] = false;
 613             bGridPBC_[ZZ] = false;
 614             break;
 615         case epbcXY:
 616             bGridPBC_[XX] = true;
 617             bGridPBC_[YY] = true;
 618             bGridPBC_[ZZ] = false;
 619             break;
 620         case epbcXYZ:
 621             bGridPBC_[XX] = true;
 622             bGridPBC_[YY] = true;
 623             bGridPBC_[ZZ] = true;
 624             break;
 625         default:
 626             // Grid searching not supported for now with screw.
 627             return false;
 628     }
 629
 630     bool   bSingleCell[DIM] = {false, false, bXY_};
 631     matrix box;
 632     copy_mat(pbc.box, box);
 633     // TODO: In principle, we could use the bounding box for periodic
 634     // dimensions as well if the bounding box is sufficiently far from the box
 635     // edges.
 636     rvec   origin, boundingBoxSize;
 637     computeBoundingBox(posCount, x, origin, boundingBoxSize);
 638     clear_rvec(gridOrigin_);
 639     for (int dd = 0; dd < DIM; ++dd)
 640     {
 641         if (!bGridPBC_[dd] && !bSingleCell[dd])
 642         {
 643             gridOrigin_[dd] = origin[dd];
 644             clear_rvec(box[dd]);
 645             box[dd][dd] = boundingBoxSize[dd];
 646         }
 647         // TODO: In case the zero vector comes from the bounding box, this does
 648         // not lead to a very efficient grid search, but that should be rare.
 649         if (box[dd][dd] <= 0.0)
 650         {
 651             GMX_ASSERT(!bGridPBC_[dd], "Periodic box vector is zero");
 652             bSingleCell[dd] = true;
 653             clear_rvec(box[dd]);
 654             box[dd][dd] = 1.0;
 655         }
 656     }
 657
 658     if (!initGridCells(box, bSingleCell, posCount))
 659     {
 660         return false;
 661     }
 662
 663     bTric_ = TRICLINIC(pbc.box);
 664     for (int dd = 0; dd < DIM; ++dd)
 665     {
 666         cellSize_[dd] = box[dd][dd] / ncelldim_[dd];
 667         if (bSingleCell[dd])
 668         {
 669             invCellSize_[dd] = 0.0;
 670         }
 671         else
 672         {
 673             invCellSize_[dd] = 1.0 / cellSize_[dd];
 674         }
 675     }
 676     if (bTric_)
 677     {
 678         cellShiftZY_ = box[ZZ][YY] * invCellSize_[YY];
 679         cellShiftZX_ = box[ZZ][XX] * invCellSize_[XX];
 680         cellShiftYX_ = box[YY][XX] * invCellSize_[XX];
 681     }
 682     return checkGridSearchEfficiency(bForce);
 683 }
 684
 685 void AnalysisNeighborhoodSearchImpl::mapPointToGridCell(const rvec x,
 686                                                         rvec       cell,
 687                                                         rvec       xout) const
 688 {
 689     rvec xtmp;
 690     rvec_sub(x, gridOrigin_, xtmp);
 691     // The reverse order is necessary for triclinic cells: shifting in Z may
 692     // modify also X and Y, and shifting in Y may modify X, so the mapping to
 693     // a rectangular grid needs to be done in this order.
 694     for (int dd = DIM - 1; dd >= 0; --dd)
 695     {
 696         real cellIndex = xtmp[dd] * invCellSize_[dd];
 697         if (bGridPBC_[dd])
 698         {
 699             const real cellCount = ncelldim_[dd];
 700             while (cellIndex < 0)
 701             {
 702                 cellIndex += cellCount;
 703                 rvec_inc(xtmp, pbc_.box[dd]);
 704             }
 705             while (cellIndex >= cellCount)
 706             {
 707                 cellIndex -= cellCount;
 708                 rvec_dec(xtmp, pbc_.box[dd]);
 709             }
 710         }
 711         cell[dd] = cellIndex;
 712     }
 713     copy_rvec(xtmp, xout);
 714 }
 715
 716 int AnalysisNeighborhoodSearchImpl::getGridCellIndex(const ivec cell) const
 717 {
 718     GMX_ASSERT(cell[XX] >= 0 && cell[XX] < ncelldim_[XX],
 719                "Grid cell X index out of range");
 720     GMX_ASSERT(cell[YY] >= 0 && cell[YY] < ncelldim_[YY],
 721                "Grid cell Y index out of range");
 722     GMX_ASSERT(cell[ZZ] >= 0 && cell[ZZ] < ncelldim_[ZZ],
 723                "Grid cell Z index out of range");
 724     return cell[XX] + cell[YY] * ncelldim_[XX]
 725            + cell[ZZ] * ncelldim_[XX] * ncelldim_[YY];
 726 }
 727
 728 void AnalysisNeighborhoodSearchImpl::addToGridCell(const rvec cell, int i)
 729 {
 730     ivec icell;
 731     for (int dd = 0; dd < DIM; ++dd)
 732     {
 733         int cellIndex = static_cast<int>(floor(cell[dd]));
 734         if (!bGridPBC_[dd])
 735         {
 736             const int cellCount = ncelldim_[dd];
 737             if (cellIndex < 0)
 738             {
 739                 cellIndex = 0;
 740             }
 741             else if (cellIndex >= cellCount)
 742             {
 743                 cellIndex = cellCount - 1;
 744             }
 745         }
 746         icell[dd] = cellIndex;
 747     }
 748     const int ci = getGridCellIndex(icell);
 749     cells_[ci].push_back(i);
 750 }
 751
 752 void AnalysisNeighborhoodSearchImpl::initCellRange(
 753         const rvec centerCell, ivec currCell, ivec upperBound, int dim) const
 754 {
 755     // TODO: Prune off cells that are completely outside the cutoff.
 756     const real range       = cutoff_ * invCellSize_[dim];
 757     real       startOffset = centerCell[dim] - range;
 758     real       endOffset   = centerCell[dim] + range;
 759     if (bTric_)
 760     {
 761         switch (dim)
 762         {
 763             case ZZ:
 764                 break;
 765             case YY:
 766                 if (currCell[ZZ] < 0)
 767                 {
 768                     startOffset += cellShiftZY_;
 769                     endOffset   += cellShiftZY_;
 770                 }
 771                 else if (currCell[ZZ] >= ncelldim_[ZZ])
 772                 {
 773                     startOffset -= cellShiftZY_;
 774                     endOffset   -= cellShiftZY_;
 775                 }
 776                 break;
 777             case XX:
 778                 if (currCell[ZZ] < 0)
 779                 {
 780                     startOffset += cellShiftZX_;
 781                     endOffset   += cellShiftZX_;
 782                 }
 783                 else if (currCell[ZZ] >= ncelldim_[ZZ])
 784                 {
 785                     startOffset -= cellShiftZX_;
 786                     endOffset   -= cellShiftZX_;
 787                 }
 788                 if (currCell[YY] < 0)
 789                 {
 790                     startOffset += cellShiftYX_;
 791                     endOffset   += cellShiftYX_;
 792                 }
 793                 else if (currCell[YY] >= ncelldim_[YY])
 794                 {
 795                     startOffset -= cellShiftYX_;
 796                     endOffset   -= cellShiftYX_;
 797                 }
 798                 break;
 799         }
 800     }
 801     // For non-periodic dimensions, clamp to the actual grid edges.
 802     if (!bGridPBC_[dim])
 803     {
 804         // If endOffset < 0 or startOffset > N, these may cause the whole
 805         // test position/grid plane/grid row to be skipped.
 806         if (startOffset < 0)
 807         {
 808             startOffset = 0;
 809         }
 810         const int cellCount = ncelldim_[dim];
 811         if (endOffset > cellCount - 1)
 812         {
 813             endOffset = cellCount - 1;
 814         }
 815     }
 816     currCell[dim]   = static_cast<int>(floor(startOffset));
 817     upperBound[dim] = static_cast<int>(floor(endOffset));
 818 }
 819
 820 bool AnalysisNeighborhoodSearchImpl::nextCell(
 821         const rvec centerCell, ivec cell, ivec upperBound) const
 822 {
 823     int dim = 0;
 824     while (dim < DIM)
 825     {
 826 next:
 827         ++cell[dim];
 828         if (cell[dim] > upperBound[dim])
 829         {
 830             ++dim;
 831             continue;
 832         }
 833         for (int d = dim - 1; d >= 0; --d)
 834         {
 835             initCellRange(centerCell, cell, upperBound, d);
 836             if (cell[d] > upperBound[d])
 837             {
 838                 dim = d + 1;
 839                 goto next;
 840             }
 841         }
 842         return true;
 843     }
 844     return false;
 845 }
 846
 847 int AnalysisNeighborhoodSearchImpl::shiftCell(const ivec cell, rvec shift) const
 848 {
 849     ivec shiftedCell;
 850     copy_ivec(cell, shiftedCell);
 851
 852     clear_rvec(shift);
 853     for (int d = 0; d < DIM; ++d)
 854     {
 855         const int cellCount = ncelldim_[d];
 856         if (bGridPBC_[d])
 857         {
 858             // A single shift may not be sufficient if the cell must be shifted
 859             // in more than one dimension, although for each individual
 860             // dimension it would be.
 861             while (shiftedCell[d] < 0)
 862             {
 863                 shiftedCell[d] += cellCount;
 864                 rvec_inc(shift, pbc_.box[d]);
 865             }
 866             while (shiftedCell[d] >= cellCount)
 867             {
 868                 shiftedCell[d] -= cellCount;
 869                 rvec_dec(shift, pbc_.box[d]);
 870             }
 871         }
 872     }
 873
 874     return getGridCellIndex(shiftedCell);
 875 }
 876
 877 void AnalysisNeighborhoodSearchImpl::init(
 878         AnalysisNeighborhood::SearchMode     mode,
 879         bool                                 bXY,
 880         const t_blocka                      *excls,
 881         const t_pbc                         *pbc,
 882         const AnalysisNeighborhoodPositions &positions)
 883 {
 884     GMX_RELEASE_ASSERT(positions.index_ == -1,
 885                        "Individual indexed positions not supported as reference");
 886     bXY_ = bXY;
 887     if (bXY_ && pbc != NULL && pbc->ePBC != epbcNONE)
 888     {
 889         if (pbc->ePBC != epbcXY && pbc->ePBC != epbcXYZ)
 890         {
 891             std::string message =
 892                 formatString("Computations in the XY plane are not supported with PBC type '%s'",
 893                              epbc_names[pbc->ePBC]);
 894             GMX_THROW(NotImplementedError(message));
 895         }
 896         if (pbc->ePBC == epbcXYZ &&
 897             (std::fabs(pbc->box[ZZ][XX]) > GMX_REAL_EPS*pbc->box[ZZ][ZZ] ||
 898              std::fabs(pbc->box[ZZ][YY]) > GMX_REAL_EPS*pbc->box[ZZ][ZZ]))
 899         {
 900             GMX_THROW(NotImplementedError("Computations in the XY plane are not supported when the last box vector is not parallel to the Z axis"));
 901         }
 902         // Use a single grid cell in Z direction.
 903         matrix box;
 904         copy_mat(pbc->box, box);
 905         clear_rvec(box[ZZ]);
 906         set_pbc(&pbc_, epbcXY, box);
 907     }
 908     else if (pbc != NULL)
 909     {
 910         pbc_ = *pbc;
 911     }
 912     else
 913     {
 914         pbc_.ePBC = epbcNONE;
 915         clear_mat(pbc_.box);
 916     }
 917     nref_ = positions.count_;
 918     if (mode == AnalysisNeighborhood::eSearchMode_Simple)
 919     {
 920         bGrid_ = false;
 921     }
 922     else if (bTryGrid_)
 923     {
 924         bGrid_ = initGrid(pbc_, positions.count_, positions.x_,
 925                           mode == AnalysisNeighborhood::eSearchMode_Grid);
 926     }
 927     refIndices_ = positions.indices_;
 928     if (bGrid_)
 929     {
 930         xrefAlloc_.resize(nref_);
 931         xref_ = as_rvec_array(xrefAlloc_.data());
 932
 933         for (int i = 0; i < nref_; ++i)
 934         {
 935             const int ii = (refIndices_ != NULL) ? refIndices_[i] : i;
 936             rvec      refcell;
 937             mapPointToGridCell(positions.x_[ii], refcell, xrefAlloc_[i]);
 938             addToGridCell(refcell, i);
 939         }
 940     }
 941     else if (refIndices_ != NULL)
 942     {
 943         xrefAlloc_.resize(nref_);
 944         xref_ = as_rvec_array(xrefAlloc_.data());
 945         for (int i = 0; i < nref_; ++i)
 946         {
 947             copy_rvec(positions.x_[refIndices_[i]], xrefAlloc_[i]);
 948         }
 949     }
 950     else
 951     {
 952         xref_ = positions.x_;
 953     }
 954     excls_           = excls;
 955     refExclusionIds_ = NULL;
 956     if (excls != NULL)
 957     {
 958         // TODO: Check that the IDs are ascending, or remove the limitation.
 959         refExclusionIds_ = positions.exclusionIds_;
 960         GMX_RELEASE_ASSERT(refExclusionIds_ != NULL,
 961                            "Exclusion IDs must be set for reference positions "
 962                            "when exclusions are enabled");
 963     }
 964 }
 965
 966 /********************************************************************
 967  * AnalysisNeighborhoodPairSearchImpl
 968  */
 969
 970 void AnalysisNeighborhoodPairSearchImpl::reset(int testIndex)
 971 {
 972     testIndex_ = testIndex;
 973     if (testIndex_ >= 0 && testIndex_ < testPosCount_)
 974     {
 975         const int index =
 976             (testIndices_ != NULL ? testIndices_[testIndex] : testIndex);
 977         if (search_.bGrid_)
 978         {
 979             search_.mapPointToGridCell(testPositions_[index], testcell_, xtest_);
 980             search_.initCellRange(testcell_, currCell_, cellBound_, ZZ);
 981             search_.initCellRange(testcell_, currCell_, cellBound_, YY);
 982             search_.initCellRange(testcell_, currCell_, cellBound_, XX);
 983         }
 984         else
 985         {
 986             copy_rvec(testPositions_[index], xtest_);
 987         }
 988         if (search_.excls_ != NULL)
 989         {
 990             const int exclIndex  = testExclusionIds_[index];
 991             if (exclIndex < search_.excls_->nr)
 992             {
 993                 const int startIndex = search_.excls_->index[exclIndex];
 994                 nexcl_ = search_.excls_->index[exclIndex + 1] - startIndex;
 995                 excl_  = &search_.excls_->a[startIndex];
 996             }
 997             else
 998             {
 999                 nexcl_ = 0;
1000                 excl_  = NULL;
1001             }
1002         }
1003     }
1004     previ_     = -1;
1005     prevr2_    = 0.0;
1006     clear_rvec(prevdx_);
1007     exclind_   = 0;
1008     prevcai_   = -1;
1009 }
1010
1011 void AnalysisNeighborhoodPairSearchImpl::nextTestPosition()
1012 {
1013     if (testIndex_ < testPosCount_)
1014     {
1015         ++testIndex_;
1016         reset(testIndex_);
1017     }
1018 }
1019
1020 bool AnalysisNeighborhoodPairSearchImpl::isExcluded(int j)
1021 {
1022     if (exclind_ < nexcl_)
1023     {
1024         const int index =
1025             (search_.refIndices_ != NULL ? search_.refIndices_[j] : j);
1026         const int refId = search_.refExclusionIds_[index];
1027         while (exclind_ < nexcl_ && excl_[exclind_] < refId)
1028         {
1029             ++exclind_;
1030         }
1031         if (exclind_ < nexcl_ && refId == excl_[exclind_])
1032         {
1033             ++exclind_;
1034             return true;
1035         }
1036     }
1037     return false;
1038 }
1039
1040 void AnalysisNeighborhoodPairSearchImpl::startSearch(
1041         const AnalysisNeighborhoodPositions &positions)
1042 {
1043     testPosCount_     = positions.count_;
1044     testPositions_    = positions.x_;
1045     testExclusionIds_ = positions.exclusionIds_;
1046     testIndices_      = positions.indices_;
1047     GMX_RELEASE_ASSERT(search_.excls_ == NULL || testExclusionIds_ != NULL,
1048                        "Exclusion IDs must be set when exclusions are enabled");
1049     if (positions.index_ < 0)
1050     {
1051         reset(0);
1052     }
1053     else
1054     {
1055         // Somewhat of a hack: setup the array such that only the last position
1056         // will be used.
1057         testPosCount_ = positions.index_ + 1;
1058         reset(positions.index_);
1059     }
1060 }
1061
1062 template <class Action>
1063 bool AnalysisNeighborhoodPairSearchImpl::searchNext(Action action)
1064 {
1065     while (testIndex_ < testPosCount_)
1066     {
1067         if (search_.bGrid_)
1068         {
1069             int cai = prevcai_ + 1;
1070
1071             do
1072             {
1073                 rvec      shift;
1074                 const int ci       = search_.shiftCell(currCell_, shift);
1075                 const int cellSize = static_cast<int>(search_.cells_[ci].size());
1076                 for (; cai < cellSize; ++cai)
1077                 {
1078                     const int i = search_.cells_[ci][cai];
1079                     if (isExcluded(i))
1080                     {
1081                         continue;
1082                     }
1083                     rvec       dx;
1084                     rvec_sub(search_.xref_[i], xtest_, dx);
1085                     rvec_sub(dx, shift, dx);
1086                     const real r2
1087                         = search_.bXY_
1088                             ? dx[XX]*dx[XX] + dx[YY]*dx[YY]
1089                             : norm2(dx);
1090                     if (r2 <= search_.cutoff2_)
1091                     {
1092                         if (action(i, r2, dx))
1093                         {
1094                             prevcai_ = cai;
1095                             previ_   = i;
1096                             prevr2_  = r2;
1097                             copy_rvec(dx, prevdx_);
1098                             return true;
1099                         }
1100                     }
1101                 }
1102                 exclind_ = 0;
1103                 cai      = 0;
1104             }
1105             while (search_.nextCell(testcell_, currCell_, cellBound_));
1106         }
1107         else
1108         {
1109             for (int i = previ_ + 1; i < search_.nref_; ++i)
1110             {
1111                 if (isExcluded(i))
1112                 {
1113                     continue;
1114                 }
1115                 rvec dx;
1116                 if (search_.pbc_.ePBC != epbcNONE)
1117                 {
1118                     pbc_dx(&search_.pbc_, search_.xref_[i], xtest_, dx);
1119                 }
1120                 else
1121                 {
1122                     rvec_sub(search_.xref_[i], xtest_, dx);
1123                 }
1124                 const real r2
1125                     = search_.bXY_
1126                         ? dx[XX]*dx[XX] + dx[YY]*dx[YY]
1127                         : norm2(dx);
1128                 if (r2 <= search_.cutoff2_)
1129                 {
1130                     if (action(i, r2, dx))
1131                     {
1132                         previ_  = i;
1133                         prevr2_ = r2;
1134                         copy_rvec(dx, prevdx_);
1135                         return true;
1136                     }
1137                 }
1138             }
1139         }
1140         nextTestPosition();
1141     }
1142     return false;
1143 }
1144
1145 void AnalysisNeighborhoodPairSearchImpl::initFoundPair(
1146         AnalysisNeighborhoodPair *pair) const
1147 {
1148     if (previ_ < 0)
1149     {
1150         *pair = AnalysisNeighborhoodPair();
1151     }
1152     else
1153     {
1154         *pair = AnalysisNeighborhoodPair(previ_, testIndex_, prevr2_, prevdx_);
1155     }
1156 }
1157
1158 }   // namespace internal
1159
1160 namespace
1161 {
1162
1163 /*! \brief
1164  * Search action to find the next neighbor.
1165  *
1166  * Used as the action for AnalysisNeighborhoodPairSearchImpl::searchNext() to
1167  * find the next neighbor.
1168  *
1169  * Simply breaks the loop on the first found neighbor.
1170  */
1171 bool withinAction(int /*i*/, real /*r2*/, const rvec /* dx */)
1172 {
1173     return true;
1174 }
1175
1176 /*! \brief
1177  * Search action find the minimum distance.
1178  *
1179  * Used as the action for AnalysisNeighborhoodPairSearchImpl::searchNext() to
1180  * find the nearest neighbor.
1181  *
1182  * With this action, AnalysisNeighborhoodPairSearchImpl::searchNext() always
1183  * returns false, and the output is put into the variables passed by pointer
1184  * into the constructor.  If no neighbors are found, the output variables are
1185  * not modified, i.e., the caller must initialize them.
1186  */
1187 class MindistAction
1188 {
1189     public:
1190         /*! \brief
1191          * Initializes the action with given output locations.
1192          *
1193          * \param[out] closestPoint Index of the closest reference location.
1194          * \param[out] minDist2     Minimum distance squared.
1195          * \param[out] dx           Shortest distance vector.
1196          *
1197          * The constructor call does not modify the pointed values, but only
1198          * stores the pointers for later use.
1199          * See the class description for additional semantics.
1200          */
1201         MindistAction(int *closestPoint, real *minDist2, rvec *dx)
1202             : closestPoint_(*closestPoint), minDist2_(*minDist2), dx_(*dx)
1203         {
1204         }
1205         //! Copies the action.
1206         MindistAction(const MindistAction &)            = default;
1207
1208         //! Processes a neighbor to find the nearest point.
1209         bool operator()(int i, real r2, const rvec dx)
1210         {
1211             if (r2 < minDist2_)
1212             {
1213                 closestPoint_ = i;
1214                 minDist2_     = r2;
1215                 copy_rvec(dx, dx_);
1216             }
1217             return false;
1218         }
1219
1220     private:
1221         int     &closestPoint_;
1222         real    &minDist2_;
1223         rvec    &dx_;
1224
1225         GMX_DISALLOW_ASSIGN(MindistAction);
1226 };
1227
1228 }   // namespace
1229
1230 /********************************************************************
1231  * AnalysisNeighborhood::Impl
1232  */
1233
1234 class AnalysisNeighborhood::Impl
1235 {
1236     public:
1237         typedef AnalysisNeighborhoodSearch::ImplPointer SearchImplPointer;
1238         typedef std::vector<SearchImplPointer> SearchList;
1239
1240         Impl()
1241             : cutoff_(0), excls_(NULL), mode_(eSearchMode_Automatic), bXY_(false)
1242         {
1243         }
1244         ~Impl()
1245         {
1246             SearchList::const_iterator i;
1247             for (i = searchList_.begin(); i != searchList_.end(); ++i)
1248             {
1249                 GMX_RELEASE_ASSERT(i->unique(),
1250                                    "Dangling AnalysisNeighborhoodSearch reference");
1251             }
1252         }
1253
1254         SearchImplPointer getSearch();
1255
1256         Mutex                   createSearchMutex_;
1257         SearchList              searchList_;
1258         real                    cutoff_;
1259         const t_blocka         *excls_;
1260         SearchMode              mode_;
1261         bool                    bXY_;
1262 };
1263
1264 AnalysisNeighborhood::Impl::SearchImplPointer
1265 AnalysisNeighborhood::Impl::getSearch()
1266 {
1267     lock_guard<Mutex> lock(createSearchMutex_);
1268     // TODO: Consider whether this needs to/can be faster, e.g., by keeping a
1269     // separate pool of unused search objects.
1270     SearchList::const_iterator i;
1271     for (i = searchList_.begin(); i != searchList_.end(); ++i)
1272     {
1273         if (i->unique())
1274         {
1275             return *i;
1276         }
1277     }
1278     SearchImplPointer search(new internal::AnalysisNeighborhoodSearchImpl(cutoff_));
1279     searchList_.push_back(search);
1280     return search;
1281 }
1282
1283 /********************************************************************
1284  * AnalysisNeighborhood
1285  */
1286
1287 AnalysisNeighborhood::AnalysisNeighborhood()
1288     : impl_(new Impl)
1289 {
1290 }
1291
1292 AnalysisNeighborhood::~AnalysisNeighborhood()
1293 {
1294 }
1295
1296 void AnalysisNeighborhood::setCutoff(real cutoff)
1297 {
1298     GMX_RELEASE_ASSERT(impl_->searchList_.empty(),
1299                        "Changing the cutoff after initSearch() not currently supported");
1300     impl_->cutoff_ = cutoff;
1301 }
1302
1303 void AnalysisNeighborhood::setXYMode(bool bXY)
1304 {
1305     impl_->bXY_ = bXY;
1306 }
1307
1308 void AnalysisNeighborhood::setTopologyExclusions(const t_blocka *excls)
1309 {
1310     GMX_RELEASE_ASSERT(impl_->searchList_.empty(),
1311                        "Changing the exclusions after initSearch() not currently supported");
1312     impl_->excls_ = excls;
1313 }
1314
1315 void AnalysisNeighborhood::setMode(SearchMode mode)
1316 {
1317     impl_->mode_ = mode;
1318 }
1319
1320 AnalysisNeighborhood::SearchMode AnalysisNeighborhood::mode() const
1321 {
1322     return impl_->mode_;
1323 }
1324
1325 AnalysisNeighborhoodSearch
1326 AnalysisNeighborhood::initSearch(const t_pbc                         *pbc,
1327                                  const AnalysisNeighborhoodPositions &positions)
1328 {
1329     Impl::SearchImplPointer search(impl_->getSearch());
1330     search->init(mode(), impl_->bXY_, impl_->excls_,
1331                  pbc, positions);
1332     return AnalysisNeighborhoodSearch(search);
1333 }
1334
1335 /********************************************************************
1336  * AnalysisNeighborhoodSearch
1337  */
1338
1339 AnalysisNeighborhoodSearch::AnalysisNeighborhoodSearch()
1340 {
1341 }
1342
1343 AnalysisNeighborhoodSearch::AnalysisNeighborhoodSearch(const ImplPointer &impl)
1344     : impl_(impl)
1345 {
1346 }
1347
1348 void AnalysisNeighborhoodSearch::reset()
1349 {
1350     impl_.reset();
1351 }
1352
1353 AnalysisNeighborhood::SearchMode AnalysisNeighborhoodSearch::mode() const
1354 {
1355     GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
1356     return (impl_->usesGridSearch()
1357             ? AnalysisNeighborhood::eSearchMode_Grid
1358             : AnalysisNeighborhood::eSearchMode_Simple);
1359 }
1360
1361 bool AnalysisNeighborhoodSearch::isWithin(
1362         const AnalysisNeighborhoodPositions &positions) const
1363 {
1364     GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
1365     internal::AnalysisNeighborhoodPairSearchImpl pairSearch(*impl_);
1366     pairSearch.startSearch(positions);
1367     return pairSearch.searchNext(&withinAction);
1368 }
1369
1370 real AnalysisNeighborhoodSearch::minimumDistance(
1371         const AnalysisNeighborhoodPositions &positions) const
1372 {
1373     GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
1374     internal::AnalysisNeighborhoodPairSearchImpl pairSearch(*impl_);
1375     pairSearch.startSearch(positions);
1376     real          minDist2     = impl_->cutoffSquared();
1377     int           closestPoint = -1;
1378     rvec          dx           = {0.0, 0.0, 0.0};
1379     MindistAction action(&closestPoint, &minDist2, &dx);
1380     (void)pairSearch.searchNext(action);
1381     return std::sqrt(minDist2);
1382 }
1383
1384 AnalysisNeighborhoodPair
1385 AnalysisNeighborhoodSearch::nearestPoint(
1386         const AnalysisNeighborhoodPositions &positions) const
1387 {
1388     GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
1389     internal::AnalysisNeighborhoodPairSearchImpl pairSearch(*impl_);
1390     pairSearch.startSearch(positions);
1391     real          minDist2     = impl_->cutoffSquared();
1392     int           closestPoint = -1;
1393     rvec          dx           = {0.0, 0.0, 0.0};
1394     MindistAction action(&closestPoint, &minDist2, &dx);
1395     (void)pairSearch.searchNext(action);
1396     return AnalysisNeighborhoodPair(closestPoint, 0, minDist2, dx);
1397 }
1398
1399 AnalysisNeighborhoodPairSearch
1400 AnalysisNeighborhoodSearch::startPairSearch(
1401         const AnalysisNeighborhoodPositions &positions) const
1402 {
1403     GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
1404     Impl::PairSearchImplPointer pairSearch(impl_->getPairSearch());
1405     pairSearch->startSearch(positions);
1406     return AnalysisNeighborhoodPairSearch(pairSearch);
1407 }
1408
1409 /********************************************************************
1410  * AnalysisNeighborhoodPairSearch
1411  */
1412
1413 AnalysisNeighborhoodPairSearch::AnalysisNeighborhoodPairSearch(
1414         const ImplPointer &impl)
1415     : impl_(impl)
1416 {
1417 }
1418
1419 bool AnalysisNeighborhoodPairSearch::findNextPair(AnalysisNeighborhoodPair *pair)
1420 {
1421     bool bFound = impl_->searchNext(&withinAction);
1422     impl_->initFoundPair(pair);
1423     return bFound;
1424 }
1425
1426 void AnalysisNeighborhoodPairSearch::skipRemainingPairsForTestPosition()
1427 {
1428     impl_->nextTestPosition();
1429 }
1430
1431 } // namespace gmx