src/gromacs/selection/nbsearch.cpp

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