2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
39 * Implements neighborhood searching for analysis (from nbsearch.h).
41 * High-level overview of the algorithm is at \ref page_analysisnbsearch.
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.
51 * \author Teemu Murtola <teemu.murtola@gmail.com>
52 * \ingroup module_selection
65 #include "gromacs/math/functions.h"
66 #include "gromacs/math/vec.h"
67 #include "gromacs/pbcutil/pbc.h"
68 #include "gromacs/utility/arrayref.h"
69 #include "gromacs/utility/classhelpers.h"
70 #include "gromacs/utility/exceptions.h"
71 #include "gromacs/utility/gmxassert.h"
72 #include "gromacs/utility/listoflists.h"
73 #include "gromacs/utility/stringutil.h"
84 * Computes the bounding box for a set of positions.
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.
91 void computeBoundingBox(int posCount, const rvec x[], rvec origin, rvec size)
94 copy_rvec(x[0], origin);
95 copy_rvec(x[0], maxBound);
96 for (int i = 1; i < posCount; ++i)
98 for (int d = 0; d < DIM; ++d)
100 if (origin[d] > x[i][d])
104 if (maxBound[d] < x[i][d])
106 maxBound[d] = x[i][d];
110 rvec_sub(maxBound, origin, size);
118 /********************************************************************
119 * Implementation class declarations
122 class AnalysisNeighborhoodSearchImpl
125 typedef AnalysisNeighborhoodPairSearch::ImplPointer PairSearchImplPointer;
126 typedef std::vector<PairSearchImplPointer> PairSearchList;
127 typedef std::vector<std::vector<int>> CellList;
129 explicit AnalysisNeighborhoodSearchImpl(real cutoff);
130 ~AnalysisNeighborhoodSearchImpl();
133 * Initializes the search with a given box and reference positions.
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.
141 void init(AnalysisNeighborhood::SearchMode mode,
143 const ListOfLists<int>* excls,
145 const AnalysisNeighborhoodPositions& positions);
146 PairSearchImplPointer getPairSearch();
148 real cutoffSquared() const { return cutoff2_; }
149 bool usesGridSearch() const { return bGrid_; }
153 * Determines a suitable grid size and sets up the cells.
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
160 * \returns `false` if grid search is not suitable.
162 bool initGridCells(const matrix box, bool bSingleCell[DIM], int posCount);
164 * Sets ua a search grid for a given box.
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.
173 bool initGrid(const t_pbc& pbc, int posCount, const rvec x[], bool bForce);
175 * Maps a point into a grid cell.
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.
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.
185 void mapPointToGridCell(const rvec x, rvec cell, rvec xout) const;
187 * Calculates linear index of a grid cell.
189 * \param[in] cell Cell indices (must be within the grid).
190 * \returns Linear index of \p cell.
192 int getGridCellIndex(const ivec cell) const;
194 * Calculates linear index of a grid cell from fractional coordinates.
196 * \param[in] cell Cell indices (must be within the grid).
197 * \returns Linear index of \p cell.
199 int getGridCellIndex(const rvec cell) const;
201 * Adds an index into a grid cell.
203 * \param[in] cell Fractional cell coordinates into which \p i should
205 * \param[in] i Index to add.
207 * \p cell should satisfy the conditions that \p mapPointToGridCell()
210 void addToGridCell(const rvec cell, int i);
212 * Initializes a cell pair loop for a dimension.
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.
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.
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
233 void initCellRange(const rvec centerCell, ivec cell, ivec upperBound, int dim) const;
235 * Computes the extent of the cutoff sphere on a particular cell edge.
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`).
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.
248 real computeCutoffExtent(RVec centerCell, const ivec cell, int dim) const;
250 * Advances cell pair loop to the next cell.
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.
257 bool nextCell(const rvec centerCell, ivec cell, ivec upperBound) const;
259 * Calculates the index and shift of a grid cell during looping.
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`.
266 int shiftCell(const ivec cell, rvec shift) const;
268 //! Whether to try grid searching.
272 //! The cutoff squared.
274 //! Whether to do searching in XY plane only.
277 //! Number of reference points for the current frame.
279 //! Reference point positions.
281 //! Reference position exclusion IDs.
282 const int* refExclusionIds_;
283 //! Reference position indices (NULL if no indices).
284 const int* refIndices_;
286 const ListOfLists<int>* excls_;
290 //! Whether grid searching is actually used for the current positions.
292 //! false if the box is rectangular.
294 //! Whether the grid is periodic in a dimension.
296 //! Array for storing in-unit-cell reference positions.
297 std::vector<RVec> xrefAlloc_;
298 //! Origin of the grid (zero for periodic dimensions).
300 //! Size of a single grid cell.
302 //! Inverse of \p cellSize_. Zero for dimensions where grid is not used.
305 * Shift in cell coordinates (for triclinic boxes) in X when crossing
306 * the Z periodic boundary.
310 * Shift in cell coordinates (for triclinic boxes) in Y when crossing
311 * the Z periodic boundary.
315 * Shift in cell coordinates (for triclinic boxes) in X when crossing
316 * the Y periodic boundary.
319 //! Number of cells along each dimension.
321 //! Data structure to hold the grid cell contents.
324 std::mutex createPairSearchMutex_;
325 PairSearchList pairSearchList_;
327 friend class AnalysisNeighborhoodPairSearchImpl;
329 GMX_DISALLOW_COPY_AND_ASSIGN(AnalysisNeighborhoodSearchImpl);
332 class AnalysisNeighborhoodPairSearchImpl
335 explicit AnalysisNeighborhoodPairSearchImpl(const AnalysisNeighborhoodSearchImpl& search) :
338 selfSearchMode_ = false;
340 testPositions_ = nullptr;
341 testExclusionIds_ = nullptr;
342 testIndices_ = nullptr;
344 clear_rvec(testcell_);
345 clear_ivec(currCell_);
346 clear_ivec(cellBound_);
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();
363 //! Clears the loop indices.
364 void reset(int testIndex);
365 //! Checks whether a reference positiong should be excluded.
366 bool isExcluded(int j);
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.
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_.
384 //! Stores test position during a pair loop.
386 //! Stores the previous returned position during a pair loop.
388 //! Stores the pair distance corresponding to previ_.
390 //! Stores the shortest distance vector corresponding to previ_.
392 //! Stores the current exclusion index during loops.
394 //! Stores the fractional test particle cell location during loops.
396 //! Stores the cell index corresponding to testcell_.
398 //! Stores the current cell during pair loops.
400 //! Stores the current loop upper bounds for each dimension during pair loops.
402 //! Stores the index within the current cell during pair loops.
405 GMX_DISALLOW_COPY_AND_ASSIGN(AnalysisNeighborhoodPairSearchImpl);
408 /********************************************************************
409 * AnalysisNeighborhoodSearchImpl
412 AnalysisNeighborhoodSearchImpl::AnalysisNeighborhoodSearchImpl(real cutoff)
418 cutoff_ = cutoff2_ = GMX_REAL_MAX;
423 cutoff2_ = gmx::square(cutoff_);
428 refExclusionIds_ = nullptr;
429 refIndices_ = nullptr;
430 std::memset(&pbc_, 0, sizeof(pbc_));
434 bGridPBC_[XX] = true;
435 bGridPBC_[YY] = true;
436 bGridPBC_[ZZ] = true;
438 clear_rvec(gridOrigin_);
439 clear_rvec(cellSize_);
440 clear_rvec(invCellSize_);
441 clear_ivec(ncelldim_);
444 AnalysisNeighborhoodSearchImpl::~AnalysisNeighborhoodSearchImpl()
446 PairSearchList::const_iterator i;
447 for (i = pairSearchList_.begin(); i != pairSearchList_.end(); ++i)
449 GMX_RELEASE_ASSERT(i->unique(), "Dangling AnalysisNeighborhoodPairSearch reference");
453 AnalysisNeighborhoodSearchImpl::PairSearchImplPointer AnalysisNeighborhoodSearchImpl::getPairSearch()
455 std::lock_guard<std::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)
466 PairSearchImplPointer pairSearch(new AnalysisNeighborhoodPairSearchImpl(*this));
467 pairSearchList_.push_back(pairSearch);
471 bool AnalysisNeighborhoodSearchImpl::initGridCells(const matrix box, bool bSingleCell[DIM], int posCount)
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;
486 for (int dd = 0; dd < DIM; ++dd)
488 const real boxSize = box[dd][dd];
489 if (boxSize < targetsize)
491 bSingleCell[dd] = true;
494 // TODO: Consider if a fallback would be possible/better.
507 if (dimCount == 0 || dimCount == prevDimCount)
511 targetsize = pow(volume * 10 / posCount, static_cast<real>(1. / dimCount));
512 prevDimCount = dimCount;
515 int totalCellCount = 1;
516 for (int dd = 0; dd < DIM; ++dd)
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)
533 totalCellCount *= cellCount;
534 ncelldim_[dd] = cellCount;
536 if (totalCellCount <= 3)
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))
545 cells_.resize(totalCellCount);
547 for (int ci = 0; ci < totalCellCount; ++ci)
554 bool AnalysisNeighborhoodSearchImpl::initGrid(const t_pbc& pbc, int posCount, const rvec x[], bool bForce)
561 // TODO: Use this again (can be useful when tuning initGridCells()),
562 // or remove throughout.
563 GMX_UNUSED_VALUE(bForce);
568 bGridPBC_[XX] = false;
569 bGridPBC_[YY] = false;
570 bGridPBC_[ZZ] = false;
573 bGridPBC_[XX] = true;
574 bGridPBC_[YY] = true;
575 bGridPBC_[ZZ] = false;
578 bGridPBC_[XX] = true;
579 bGridPBC_[YY] = true;
580 bGridPBC_[ZZ] = true;
583 // Grid searching not supported for now with screw.
587 bool bSingleCell[DIM] = { false, false, bXY_ };
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
593 rvec origin, boundingBoxSize;
594 computeBoundingBox(posCount, x, origin, boundingBoxSize);
595 clear_rvec(gridOrigin_);
596 for (int dd = 0; dd < DIM; ++dd)
598 if (!bGridPBC_[dd] && !bSingleCell[dd])
600 gridOrigin_[dd] = origin[dd];
602 box[dd][dd] = boundingBoxSize[dd];
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)
608 GMX_ASSERT(!bGridPBC_[dd], "Periodic box vector is zero");
609 bSingleCell[dd] = true;
615 if (!initGridCells(box, bSingleCell, posCount))
620 bTric_ = TRICLINIC(pbc.box);
621 for (int dd = 0; dd < DIM; ++dd)
623 cellSize_[dd] = box[dd][dd] / ncelldim_[dd];
626 invCellSize_[dd] = 0.0;
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])
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
645 cellShiftZY_ = box[ZZ][YY] * invCellSize_[YY];
646 cellShiftZX_ = box[ZZ][XX] * invCellSize_[XX];
647 cellShiftYX_ = box[YY][XX] * invCellSize_[XX];
652 void AnalysisNeighborhoodSearchImpl::mapPointToGridCell(const rvec x, rvec cell, rvec xout) const
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)
661 real cellIndex = xtmp[dd] * invCellSize_[dd];
664 const real cellCount = ncelldim_[dd];
665 while (cellIndex < 0)
667 cellIndex += cellCount;
668 rvec_inc(xtmp, pbc_.box[dd]);
670 while (cellIndex >= cellCount)
672 cellIndex -= cellCount;
673 rvec_dec(xtmp, pbc_.box[dd]);
676 cell[dd] = cellIndex;
678 copy_rvec(xtmp, xout);
681 int AnalysisNeighborhoodSearchImpl::getGridCellIndex(const ivec cell) const
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];
689 int AnalysisNeighborhoodSearchImpl::getGridCellIndex(const rvec cell) const
692 for (int dd = 0; dd < DIM; ++dd)
694 int cellIndex = static_cast<int>(floor(cell[dd]));
697 const int cellCount = ncelldim_[dd];
702 else if (cellIndex >= cellCount)
704 cellIndex = cellCount - 1;
707 icell[dd] = cellIndex;
709 return getGridCellIndex(icell);
712 void AnalysisNeighborhoodSearchImpl::addToGridCell(const rvec cell, int i)
714 const int ci = getGridCellIndex(cell);
715 cells_[ci].push_back(i);
718 void AnalysisNeighborhoodSearchImpl::initCellRange(const rvec centerCell, ivec currCell, ivec upperBound, int dim) const
720 RVec shiftedCenter(centerCell);
721 // Shift the center to the cell coordinates of currCell, so that
722 // computeCutoffExtent() can assume simple rectangular grid.
727 if (currCell[ZZ] < 0)
729 shiftedCenter[XX] += cellShiftZX_;
731 else if (currCell[ZZ] >= ncelldim_[ZZ])
733 shiftedCenter[XX] -= cellShiftZX_;
735 if (currCell[YY] < 0)
737 shiftedCenter[XX] += cellShiftYX_;
739 else if (currCell[YY] >= ncelldim_[YY])
741 shiftedCenter[XX] -= cellShiftYX_;
744 if (dim == XX || dim == YY)
746 if (currCell[ZZ] < 0)
748 shiftedCenter[YY] += cellShiftZY_;
750 else if (currCell[ZZ] >= ncelldim_[ZZ])
752 shiftedCenter[YY] -= cellShiftZY_;
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.
762 // If endOffset < 0 or startOffset > N, these may cause the whole
763 // test position/grid plane/grid row to be skipped.
768 const int cellCount = ncelldim_[dim];
769 if (endOffset > cellCount - 1)
771 endOffset = cellCount - 1;
774 currCell[dim] = static_cast<int>(floor(startOffset));
775 upperBound[dim] = static_cast<int>(floor(endOffset));
778 real AnalysisNeighborhoodSearchImpl::computeCutoffExtent(const RVec centerCell, const ivec cell, int dim) const
786 for (int d = dim + 1; d < DIM; ++d)
788 real dimDist = cell[d] - centerCell[d];
793 else if (dimDist <= 0)
797 dist2 += dimDist * dimDist * cellSize_[d] * cellSize_[d];
799 if (dist2 >= cutoff2_)
803 return std::sqrt(cutoff2_ - dist2);
806 bool AnalysisNeighborhoodSearchImpl::nextCell(const rvec centerCell, ivec cell, ivec upperBound) const
813 if (cell[dim] > upperBound[dim])
818 for (int d = dim - 1; d >= 0; --d)
820 initCellRange(centerCell, cell, upperBound, d);
821 if (cell[d] > upperBound[d])
832 int AnalysisNeighborhoodSearchImpl::shiftCell(const ivec cell, rvec shift) const
835 copy_ivec(cell, shiftedCell);
838 for (int d = 0; d < DIM; ++d)
840 const int cellCount = ncelldim_[d];
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)
848 shiftedCell[d] += cellCount;
849 rvec_inc(shift, pbc_.box[d]);
851 while (shiftedCell[d] >= cellCount)
853 shiftedCell[d] -= cellCount;
854 rvec_dec(shift, pbc_.box[d]);
859 return getGridCellIndex(shiftedCell);
862 void AnalysisNeighborhoodSearchImpl::init(AnalysisNeighborhood::SearchMode mode,
864 const ListOfLists<int>* excls,
866 const AnalysisNeighborhoodPositions& positions)
868 GMX_RELEASE_ASSERT(positions.index_ == -1,
869 "Individual indexed positions not supported as reference");
871 if (bXY_ && pbc != nullptr && pbc->pbcType != PbcType::No)
873 if (pbc->pbcType != PbcType::XY && pbc->pbcType != PbcType::Xyz)
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));
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]))
885 NotImplementedError("Computations in the XY plane are not supported when the "
886 "last box vector is not parallel to the Z axis"));
888 // Use a single grid cell in Z direction.
890 copy_mat(pbc->box, box);
892 set_pbc(&pbc_, PbcType::XY, box);
894 else if (pbc != nullptr)
900 pbc_.pbcType = PbcType::No;
903 nref_ = positions.count_;
904 if (mode == AnalysisNeighborhood::eSearchMode_Simple)
911 pbc_, positions.count_, positions.x_, mode == AnalysisNeighborhood::eSearchMode_Grid);
913 refIndices_ = positions.indices_;
916 xrefAlloc_.resize(nref_);
917 xref_ = as_rvec_array(xrefAlloc_.data());
919 for (int i = 0; i < nref_; ++i)
921 const int ii = (refIndices_ != nullptr) ? refIndices_[i] : i;
923 mapPointToGridCell(positions.x_[ii], refcell, xrefAlloc_[i]);
924 addToGridCell(refcell, i);
927 else if (refIndices_ != nullptr)
929 xrefAlloc_.resize(nref_);
930 xref_ = as_rvec_array(xrefAlloc_.data());
931 for (int i = 0; i < nref_; ++i)
933 copy_rvec(positions.x_[refIndices_[i]], xrefAlloc_[i]);
938 xref_ = positions.x_;
941 refExclusionIds_ = nullptr;
942 if (excls != nullptr)
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");
952 /********************************************************************
953 * AnalysisNeighborhoodPairSearchImpl
956 void AnalysisNeighborhoodPairSearchImpl::reset(int testIndex)
958 testIndex_ = testIndex;
965 if (testIndex_ >= 0 && testIndex_ < testPosCount_)
967 const int index = (testIndices_ != nullptr ? testIndices_[testIndex] : testIndex);
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);
976 testCellIndex_ = search_.getGridCellIndex(testcell_);
981 copy_rvec(testPositions_[index], xtest_);
987 if (search_.excls_ != nullptr)
989 const int exclIndex = testExclusionIds_[index];
990 if (exclIndex < search_.excls_->ssize())
992 excl_ = (*search_.excls_)[exclIndex];
996 excl_ = ArrayRef<const int>();
1002 void AnalysisNeighborhoodPairSearchImpl::nextTestPosition()
1004 if (testIndex_ < testPosCount_)
1011 bool AnalysisNeighborhoodPairSearchImpl::isExcluded(int j)
1013 const int nexcl = excl_.ssize();
1014 if (exclind_ < nexcl)
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)
1022 if (exclind_ < nexcl && refId == excl_[exclind_])
1031 void AnalysisNeighborhoodPairSearchImpl::startSearch(const AnalysisNeighborhoodPositions& positions)
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)
1046 // Somewhat of a hack: setup the array such that only the last position
1048 testPosCount_ = positions.index_ + 1;
1049 reset(positions.index_);
1053 void AnalysisNeighborhoodPairSearchImpl::startSelfSearch()
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");
1065 template<class Action>
1066 bool AnalysisNeighborhoodPairSearchImpl::searchNext(Action action)
1068 while (testIndex_ < testPosCount_)
1072 int cai = prevcai_ + 1;
1077 const int ci = search_.shiftCell(currCell_, shift);
1078 if (selfSearchMode_ && ci > testCellIndex_)
1082 const int cellSize = ssize(search_.cells_[ci]);
1083 for (; cai < cellSize; ++cai)
1085 const int i = search_.cells_[ci][cai];
1086 if (selfSearchMode_ && ci == testCellIndex_ && i >= testIndex_)
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_)
1100 if (action(i, r2, dx))
1105 copy_rvec(dx, prevdx_);
1112 } while (search_.nextCell(testcell_, currCell_, cellBound_));
1116 for (int i = previ_ + 1; i < search_.nref_; ++i)
1123 if (search_.pbc_.pbcType != PbcType::No)
1125 pbc_dx(&search_.pbc_, search_.xref_[i], xtest_, dx);
1129 rvec_sub(search_.xref_[i], xtest_, dx);
1131 const real r2 = search_.bXY_ ? dx[XX] * dx[XX] + dx[YY] * dx[YY] : norm2(dx);
1132 if (r2 <= search_.cutoff2_)
1134 if (action(i, r2, dx))
1138 copy_rvec(dx, prevdx_);
1149 void AnalysisNeighborhoodPairSearchImpl::initFoundPair(AnalysisNeighborhoodPair* pair) const
1153 *pair = AnalysisNeighborhoodPair();
1157 *pair = AnalysisNeighborhoodPair(previ_, testIndex_, prevr2_, prevdx_);
1161 } // namespace internal
1167 * Search action to find the next neighbor.
1169 * Used as the action for AnalysisNeighborhoodPairSearchImpl::searchNext() to
1170 * find the next neighbor.
1172 * Simply breaks the loop on the first found neighbor.
1174 bool withinAction(int /*i*/, real /*r2*/, const rvec /* dx */)
1180 * Search action find the minimum distance.
1182 * Used as the action for AnalysisNeighborhoodPairSearchImpl::searchNext() to
1183 * find the nearest neighbor.
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.
1194 * Initializes the action with given output locations.
1196 * \param[out] closestPoint Index of the closest reference location.
1197 * \param[out] minDist2 Minimum distance squared.
1198 * \param[out] dx Shortest distance vector.
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.
1204 MindistAction(int* closestPoint, real* minDist2, rvec* dx) // NOLINT(readability-non-const-parameter)
1206 closestPoint_(*closestPoint), minDist2_(*minDist2), dx_(*dx)
1209 //! Copies the action.
1210 MindistAction(const MindistAction&) = default;
1212 //! Processes a neighbor to find the nearest point.
1213 bool operator()(int i, real r2, const rvec dx)
1229 GMX_DISALLOW_ASSIGN(MindistAction);
1234 /********************************************************************
1235 * AnalysisNeighborhood::Impl
1238 class AnalysisNeighborhood::Impl
1241 typedef AnalysisNeighborhoodSearch::ImplPointer SearchImplPointer;
1242 typedef std::vector<SearchImplPointer> SearchList;
1244 Impl() : cutoff_(0), excls_(nullptr), mode_(eSearchMode_Automatic), bXY_(false) {}
1247 SearchList::const_iterator i;
1248 for (i = searchList_.begin(); i != searchList_.end(); ++i)
1250 GMX_RELEASE_ASSERT(i->unique(), "Dangling AnalysisNeighborhoodSearch reference");
1254 SearchImplPointer getSearch();
1256 std::mutex createSearchMutex_;
1257 SearchList searchList_;
1259 const ListOfLists<int>* excls_;
1264 AnalysisNeighborhood::Impl::SearchImplPointer AnalysisNeighborhood::Impl::getSearch()
1266 std::lock_guard<std::mutex> lock(createSearchMutex_);
1267 // TODO: Consider whether this needs to/can be faster, e.g., by keeping a
1268 // separate pool of unused search objects.
1269 SearchList::const_iterator i;
1270 for (i = searchList_.begin(); i != searchList_.end(); ++i)
1277 SearchImplPointer search(new internal::AnalysisNeighborhoodSearchImpl(cutoff_));
1278 searchList_.push_back(search);
1282 /********************************************************************
1283 * AnalysisNeighborhood
1286 AnalysisNeighborhood::AnalysisNeighborhood() : impl_(new Impl) {}
1288 AnalysisNeighborhood::~AnalysisNeighborhood() {}
1290 void AnalysisNeighborhood::setCutoff(real cutoff)
1292 GMX_RELEASE_ASSERT(impl_->searchList_.empty(),
1293 "Changing the cutoff after initSearch() not currently supported");
1294 impl_->cutoff_ = cutoff;
1297 void AnalysisNeighborhood::setXYMode(bool bXY)
1302 void AnalysisNeighborhood::setTopologyExclusions(const ListOfLists<int>* excls)
1304 GMX_RELEASE_ASSERT(impl_->searchList_.empty(),
1305 "Changing the exclusions after initSearch() not currently supported");
1306 impl_->excls_ = excls;
1309 void AnalysisNeighborhood::setMode(SearchMode mode)
1311 impl_->mode_ = mode;
1314 AnalysisNeighborhood::SearchMode AnalysisNeighborhood::mode() const
1316 return impl_->mode_;
1319 AnalysisNeighborhoodSearch AnalysisNeighborhood::initSearch(const t_pbc* pbc,
1320 const AnalysisNeighborhoodPositions& positions)
1322 Impl::SearchImplPointer search(impl_->getSearch());
1323 search->init(mode(), impl_->bXY_, impl_->excls_, pbc, positions);
1324 return AnalysisNeighborhoodSearch(search);
1327 /********************************************************************
1328 * AnalysisNeighborhoodSearch
1331 AnalysisNeighborhoodSearch::AnalysisNeighborhoodSearch() {}
1333 AnalysisNeighborhoodSearch::AnalysisNeighborhoodSearch(const ImplPointer& impl) : impl_(impl) {}
1335 void AnalysisNeighborhoodSearch::reset()
1340 AnalysisNeighborhood::SearchMode AnalysisNeighborhoodSearch::mode() const
1342 GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
1343 return (impl_->usesGridSearch() ? AnalysisNeighborhood::eSearchMode_Grid
1344 : AnalysisNeighborhood::eSearchMode_Simple);
1347 bool AnalysisNeighborhoodSearch::isWithin(const AnalysisNeighborhoodPositions& positions) const
1349 GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
1350 internal::AnalysisNeighborhoodPairSearchImpl pairSearch(*impl_);
1351 pairSearch.startSearch(positions);
1352 return pairSearch.searchNext(&withinAction);
1355 real AnalysisNeighborhoodSearch::minimumDistance(const AnalysisNeighborhoodPositions& positions) const
1357 GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
1358 internal::AnalysisNeighborhoodPairSearchImpl pairSearch(*impl_);
1359 pairSearch.startSearch(positions);
1360 real minDist2 = impl_->cutoffSquared();
1361 int closestPoint = -1;
1362 rvec dx = { 0.0, 0.0, 0.0 };
1363 MindistAction action(&closestPoint, &minDist2, &dx);
1364 (void)pairSearch.searchNext(action);
1365 return std::sqrt(minDist2);
1368 AnalysisNeighborhoodPair AnalysisNeighborhoodSearch::nearestPoint(const AnalysisNeighborhoodPositions& positions) const
1370 GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
1371 internal::AnalysisNeighborhoodPairSearchImpl pairSearch(*impl_);
1372 pairSearch.startSearch(positions);
1373 real minDist2 = impl_->cutoffSquared();
1374 int closestPoint = -1;
1375 rvec dx = { 0.0, 0.0, 0.0 };
1376 MindistAction action(&closestPoint, &minDist2, &dx);
1377 (void)pairSearch.searchNext(action);
1378 return AnalysisNeighborhoodPair(closestPoint, 0, minDist2, dx);
1381 AnalysisNeighborhoodPairSearch AnalysisNeighborhoodSearch::startSelfPairSearch() const
1383 GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
1384 Impl::PairSearchImplPointer pairSearch(impl_->getPairSearch());
1385 pairSearch->startSelfSearch();
1386 return AnalysisNeighborhoodPairSearch(pairSearch);
1389 AnalysisNeighborhoodPairSearch
1390 AnalysisNeighborhoodSearch::startPairSearch(const AnalysisNeighborhoodPositions& positions) const
1392 GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
1393 Impl::PairSearchImplPointer pairSearch(impl_->getPairSearch());
1394 pairSearch->startSearch(positions);
1395 return AnalysisNeighborhoodPairSearch(pairSearch);
1398 /********************************************************************
1399 * AnalysisNeighborhoodPairSearch
1402 AnalysisNeighborhoodPairSearch::AnalysisNeighborhoodPairSearch(const ImplPointer& impl) :
1407 bool AnalysisNeighborhoodPairSearch::findNextPair(AnalysisNeighborhoodPair* pair)
1409 bool bFound = impl_->searchNext(&withinAction);
1410 impl_->initFoundPair(pair);
1414 void AnalysisNeighborhoodPairSearch::skipRemainingPairsForTestPosition()
1416 impl_->nextTestPosition();