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