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