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