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