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