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