Tidy: modernize-use-nullptr
[gromacs.git] / src / gromacs / selection / nbsearch.cpp
blob9dcfa19d203ead8719ea373d4f2671200fe2a35c
1 /*
2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2009,2010,2011,2012,2013,2014,2015,2016,2017, 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.
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.
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.
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.
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.
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.
35 /*! \internal \file
36 * \brief
37 * Implements neighborhood searching for analysis (from nbsearch.h).
39 * High-level overview of the algorithm is at \ref page_analysisnbsearch.
41 * \todo
42 * The grid implementation could still be optimized in several different ways:
43 * - Pruning grid cells from the search list if they are completely outside
44 * the sphere that is being considered.
45 * - A better heuristic could be added for falling back to simple loops for a
46 * small number of reference particles.
47 * - A better heuristic for selecting the grid size.
48 * - A multi-level grid implementation could be used to be able to use small
49 * grids for short cutoffs with very inhomogeneous particle distributions
50 * without a memory cost.
52 * \author Teemu Murtola <teemu.murtola@gmail.com>
53 * \ingroup module_selection
55 #include "gmxpre.h"
57 #include "nbsearch.h"
59 #include <cmath>
60 #include <cstring>
62 #include <algorithm>
63 #include <vector>
65 #include "gromacs/math/functions.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/mutex.h"
74 #include "gromacs/utility/stringutil.h"
76 namespace gmx
79 namespace
82 /*! \brief
83 * Computes the bounding box for a set of positions.
85 * \param[in] posCount Number of positions in \p x.
86 * \param[in] x Positions to compute the bounding box for.
87 * \param[out] origin Origin of the bounding box.
88 * \param[out] size Size of the bounding box.
90 void computeBoundingBox(int posCount, const rvec x[], rvec origin, rvec size)
92 rvec maxBound;
93 copy_rvec(x[0], origin);
94 copy_rvec(x[0], maxBound);
95 for (int i = 1; i < posCount; ++i)
97 for (int d = 0; d < DIM; ++d)
99 if (origin[d] > x[i][d])
101 origin[d] = x[i][d];
103 if (maxBound[d] < x[i][d])
105 maxBound[d] = x[i][d];
109 rvec_sub(maxBound, origin, size);
112 } // namespace
114 namespace internal
117 /********************************************************************
118 * Implementation class declarations
121 class AnalysisNeighborhoodSearchImpl
123 public:
124 typedef AnalysisNeighborhoodPairSearch::ImplPointer
125 PairSearchImplPointer;
126 typedef std::vector<PairSearchImplPointer> PairSearchList;
127 typedef std::vector<std::vector<int> > CellList;
129 explicit AnalysisNeighborhoodSearchImpl(real cutoff);
130 ~AnalysisNeighborhoodSearchImpl();
132 /*! \brief
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,
142 bool bXY,
143 const t_blocka *excls,
144 const t_pbc *pbc,
145 const AnalysisNeighborhoodPositions &positions);
146 PairSearchImplPointer getPairSearch();
148 real cutoffSquared() const { return cutoff2_; }
149 bool usesGridSearch() const { return bGrid_; }
151 private:
152 /*! \brief
153 * Checks the efficiency and possibility of doing grid-based searching.
155 * \param[in] bForce If `true`, grid search will be forced if possible.
156 * \returns `false` if grid search is not suitable.
158 bool checkGridSearchEfficiency(bool bForce);
159 /*! \brief
160 * Determines a suitable grid size and sets up the cells.
162 * \param[in] box Box vectors (should not have zero vectors).
163 * \param[in] bSingleCell If `true`, the corresponding dimension will
164 * be forced to use a single cell.
165 * \param[in] posCount Number of positions that will be put on the
166 * grid.
167 * \returns `false` if grid search is not suitable.
169 bool initGridCells(const matrix box, bool bSingleCell[DIM],
170 int posCount);
171 /*! \brief
172 * Sets ua a search grid for a given box.
174 * \param[in] pbc Information about the box.
175 * \param[in] posCount Number of positions in \p x.
176 * \param[in] x Reference positions that will be put on the grid.
177 * \param[in] bForce If `true`, grid searching will be used if at all
178 * possible, even if a simple search might give better performance.
179 * \returns `false` if grid search is not suitable.
181 bool initGrid(const t_pbc &pbc, int posCount, const rvec x[], bool bForce);
182 /*! \brief
183 * Maps a point into a grid cell.
185 * \param[in] x Point to map.
186 * \param[out] cell Fractional cell coordinates of \p x on the grid.
187 * \param[out] xout Coordinates to use.
189 * \p xout will be within the rectangular unit cell in dimensions where
190 * the grid is periodic. For other dimensions, both \p xout and
191 * \p cell can be outside the grid/unit cell.
193 void mapPointToGridCell(const rvec x, rvec cell, rvec xout) const;
194 /*! \brief
195 * Calculates linear index of a grid cell.
197 * \param[in] cell Cell indices (must be within the grid).
198 * \returns Linear index of \p cell.
200 int getGridCellIndex(const ivec cell) const;
201 /*! \brief
202 * Adds an index into a grid cell.
204 * \param[in] cell Fractional cell coordinates into which \p i should
205 * be added.
206 * \param[in] i Index to add.
208 * \p cell should satisfy the conditions that \p mapPointToGridCell()
209 * produces.
211 void addToGridCell(const rvec cell, int i);
212 /*! \brief
213 * Initializes a cell pair loop for a dimension.
215 * \param[in] centerCell Fractional cell coordiates of the particle
216 * for which pairs are being searched.
217 * \param[in,out] cell Current/initial cell to loop over.
218 * \param[in,out] upperBound Last cell to loop over.
219 * \param[in] dim Dimension to initialize in this call.
221 * Initializes `cell[dim]` and `upperBound[dim]` for looping over
222 * neighbors of a particle at position given by \p centerCell.
223 * If 'dim != ZZ`, `cell[d]` (`d > dim`) set the plane/row of cells
224 * for which the loop is initialized. The loop should then go from
225 * `cell[dim]` until `upperBound[dim]`, inclusive.
226 * `cell[d]` with `d < dim` or `upperBound[d]` with `d != dim` are not
227 * modified by this function.
229 * `cell` and `upperBound` may be outside the grid for periodic
230 * dimensions and need to be shifted separately: to simplify the
231 * looping, the range is always (roughly) symmetric around the value in
232 * `centerCell`.
234 void initCellRange(const rvec centerCell, ivec cell,
235 ivec upperBound, int dim) const;
236 /*! \brief
237 * Advances cell pair loop to the next cell.
239 * \param[in] centerCell Fractional cell coordiates of the particle
240 * for which pairs are being searched.
241 * \param[in,out] cell Current (in)/next (out) cell in the loop.
242 * \param[in,out] upperBound Last cell in the loop for each dimension.
244 bool nextCell(const rvec centerCell, ivec cell, ivec upperBound) const;
245 /*! \brief
246 * Calculates the index and shift of a grid cell during looping.
248 * \param[in] cell Unshifted cell index.
249 * \param[out] shift Shift to apply to get the periodic distance
250 * for distances between the cells.
251 * \returns Grid cell index corresponding to `cell`.
253 int shiftCell(const ivec cell, rvec shift) const;
255 //! Whether to try grid searching.
256 bool bTryGrid_;
257 //! The cutoff.
258 real cutoff_;
259 //! The cutoff squared.
260 real cutoff2_;
261 //! Whether to do searching in XY plane only.
262 bool bXY_;
264 //! Number of reference points for the current frame.
265 int nref_;
266 //! Reference point positions.
267 const rvec *xref_;
268 //! Reference position exclusion IDs.
269 const int *refExclusionIds_;
270 //! Reference position indices (NULL if no indices).
271 const int *refIndices_;
272 //! Exclusions.
273 const t_blocka *excls_;
274 //! PBC data.
275 t_pbc pbc_;
277 //! Whether grid searching is actually used for the current positions.
278 bool bGrid_;
279 //! false if the box is rectangular.
280 bool bTric_;
281 //! Whether the grid is periodic in a dimension.
282 bool bGridPBC_[DIM];
283 //! Array for storing in-unit-cell reference positions.
284 std::vector<RVec> xrefAlloc_;
285 //! Origin of the grid (zero for periodic dimensions).
286 rvec gridOrigin_;
287 //! Size of a single grid cell.
288 rvec cellSize_;
289 //! Inverse of \p cellSize_. Zero for dimensions where grid is not used.
290 rvec invCellSize_;
291 /*! \brief
292 * Shift in cell coordinates (for triclinic boxes) in X when crossing
293 * the Z periodic boundary.
295 real cellShiftZX_;
296 /*! \brief
297 * Shift in cell coordinates (for triclinic boxes) in Y when crossing
298 * the Z periodic boundary.
300 real cellShiftZY_;
301 /*! \brief
302 * Shift in cell coordinates (for triclinic boxes) in X when crossing
303 * the Y periodic boundary.
305 real cellShiftYX_;
306 //! Number of cells along each dimension.
307 ivec ncelldim_;
308 //! Data structure to hold the grid cell contents.
309 CellList cells_;
311 Mutex createPairSearchMutex_;
312 PairSearchList pairSearchList_;
314 friend class AnalysisNeighborhoodPairSearchImpl;
316 GMX_DISALLOW_COPY_AND_ASSIGN(AnalysisNeighborhoodSearchImpl);
319 class AnalysisNeighborhoodPairSearchImpl
321 public:
322 explicit AnalysisNeighborhoodPairSearchImpl(const AnalysisNeighborhoodSearchImpl &search)
323 : search_(search)
325 testPosCount_ = 0;
326 testPositions_ = nullptr;
327 testExclusionIds_ = nullptr;
328 testIndices_ = nullptr;
329 nexcl_ = 0;
330 excl_ = nullptr;
331 clear_rvec(xtest_);
332 clear_rvec(testcell_);
333 clear_ivec(currCell_);
334 clear_ivec(cellBound_);
335 reset(-1);
338 //! Initializes a search to find reference positions neighboring \p x.
339 void startSearch(const AnalysisNeighborhoodPositions &positions);
340 //! Searches for the next neighbor.
341 template <class Action>
342 bool searchNext(Action action);
343 //! Initializes a pair representing the pair found by searchNext().
344 void initFoundPair(AnalysisNeighborhoodPair *pair) const;
345 //! Advances to the next test position, skipping any remaining pairs.
346 void nextTestPosition();
348 private:
349 //! Clears the loop indices.
350 void reset(int testIndex);
351 //! Checks whether a reference positiong should be excluded.
352 bool isExcluded(int j);
354 //! Parent search object.
355 const AnalysisNeighborhoodSearchImpl &search_;
356 //! Number of test positions.
357 int testPosCount_;
358 //! Reference to the test positions.
359 const rvec *testPositions_;
360 //! Reference to the test exclusion indices.
361 const int *testExclusionIds_;
362 //! Reference to the test position indices.
363 const int *testIndices_;
364 //! Number of excluded reference positions for current test particle.
365 int nexcl_;
366 //! Exclusions for current test particle.
367 const int *excl_;
368 //! Index of the currently active test position in \p testPositions_.
369 int testIndex_;
370 //! Stores test position during a pair loop.
371 rvec xtest_;
372 //! Stores the previous returned position during a pair loop.
373 int previ_;
374 //! Stores the pair distance corresponding to previ_;
375 real prevr2_;
376 //! Stores the shortest distance vector corresponding to previ_;
377 rvec prevdx_;
378 //! Stores the current exclusion index during loops.
379 int exclind_;
380 //! Stores the fractional test particle cell location during loops.
381 rvec testcell_;
382 //! Stores the current cell during pair loops.
383 ivec currCell_;
384 //! Stores the current loop upper bounds for each dimension during pair loops.
385 ivec cellBound_;
386 //! Stores the index within the current cell during pair loops.
387 int prevcai_;
389 GMX_DISALLOW_COPY_AND_ASSIGN(AnalysisNeighborhoodPairSearchImpl);
392 /********************************************************************
393 * AnalysisNeighborhoodSearchImpl
396 AnalysisNeighborhoodSearchImpl::AnalysisNeighborhoodSearchImpl(real cutoff)
398 bTryGrid_ = true;
399 cutoff_ = cutoff;
400 if (cutoff_ <= 0)
402 cutoff_ = cutoff2_ = GMX_REAL_MAX;
403 bTryGrid_ = false;
405 else
407 cutoff2_ = gmx::square(cutoff_);
409 bXY_ = false;
410 nref_ = 0;
411 xref_ = nullptr;
412 refExclusionIds_ = nullptr;
413 refIndices_ = nullptr;
414 std::memset(&pbc_, 0, sizeof(pbc_));
416 bGrid_ = false;
417 bTric_ = false;
418 bGridPBC_[XX] = true;
419 bGridPBC_[YY] = true;
420 bGridPBC_[ZZ] = true;
422 clear_rvec(gridOrigin_);
423 clear_rvec(cellSize_);
424 clear_rvec(invCellSize_);
425 clear_ivec(ncelldim_);
428 AnalysisNeighborhoodSearchImpl::~AnalysisNeighborhoodSearchImpl()
430 PairSearchList::const_iterator i;
431 for (i = pairSearchList_.begin(); i != pairSearchList_.end(); ++i)
433 GMX_RELEASE_ASSERT(i->unique(),
434 "Dangling AnalysisNeighborhoodPairSearch reference");
438 AnalysisNeighborhoodSearchImpl::PairSearchImplPointer
439 AnalysisNeighborhoodSearchImpl::getPairSearch()
441 lock_guard<Mutex> lock(createPairSearchMutex_);
442 // TODO: Consider whether this needs to/can be faster, e.g., by keeping a
443 // separate pool of unused search objects.
444 PairSearchList::const_iterator i;
445 for (i = pairSearchList_.begin(); i != pairSearchList_.end(); ++i)
447 if (i->unique())
449 return *i;
452 PairSearchImplPointer pairSearch(new AnalysisNeighborhoodPairSearchImpl(*this));
453 pairSearchList_.push_back(pairSearch);
454 return pairSearch;
457 bool AnalysisNeighborhoodSearchImpl::checkGridSearchEfficiency(bool bForce)
459 // Find the extent of the sphere in cells.
460 ivec range;
461 for (int dd = 0; dd < DIM; ++dd)
463 range[dd] = static_cast<int>(ceil(cutoff_ * invCellSize_[dd]));
466 // Calculate the fraction of cell pairs that need to be searched,
467 // and check that the cutoff is not too large for periodic dimensions.
468 real coveredCells = 1.0;
469 for (int dd = 0; dd < DIM; ++dd)
471 const int cellCount = ncelldim_[dd];
472 const int coveredCount = 2 * range[dd] + 1;
473 if (bGridPBC_[dd])
475 if (coveredCount > cellCount)
477 // Cutoff is too close to half the box size for grid searching
478 // (it is not possible to find a single shift for every pair of
479 // grid cells).
480 return false;
482 coveredCells *= coveredCount;
484 else
486 if (range[dd] >= cellCount - 1)
488 range[dd] = cellCount - 1;
489 coveredCells *= cellCount;
491 else if (coveredCount > cellCount)
493 // The sum of range+1, range+2, ..., range+N/2, ... range+1.
494 coveredCells *= range[dd] +
495 static_cast<real>((cellCount + 1)/2 * (cellCount/2 + 1)) / cellCount;
497 else
499 // The sum of range+1, ..., 2*range+1, ..., 2*range+1, ... range+1.
500 coveredCells *= coveredCount -
501 static_cast<real>(range[dd] * (range[dd] + 1)) / cellCount;
505 // Magic constant that would need tuning for optimal performance:
506 // Don't do grid searching if nearly all cell pairs would anyways need to
507 // be looped through.
508 const int totalCellCount = ncelldim_[XX] * ncelldim_[YY] * ncelldim_[ZZ];
509 if (!bForce && coveredCells >= 0.5 * totalCellCount)
511 return false;
513 return true;
516 bool AnalysisNeighborhoodSearchImpl::initGridCells(
517 const matrix box, bool bSingleCell[DIM], int posCount)
519 // Determine the size of cubes where there are on average 10 positions.
520 // The loop takes care of cases where some of the box edges are shorter
521 // than the the desired cube size; in such cases, a single grid cell is
522 // used in these dimensions, and the cube size is determined only from the
523 // larger box vectors. Such boxes should be rare, but the bounding box
524 // approach can result in very flat boxes with certain types of selections
525 // (e.g., for interfacial systems or for small number of atoms).
526 real targetsize = 0.0;
527 int prevDimCount = 4;
528 while (true)
530 real volume = 1.0;
531 int dimCount = 3;
532 for (int dd = 0; dd < DIM; ++dd)
534 const real boxSize = box[dd][dd];
535 if (boxSize < targetsize)
537 bSingleCell[dd] = true;
538 if (bGridPBC_[dd])
540 return false;
543 if (bSingleCell[dd])
545 --dimCount;
547 else
549 volume *= boxSize;
552 if (dimCount == 0 || dimCount == prevDimCount)
554 break;
556 targetsize = pow(volume * 10 / posCount, static_cast<real>(1./dimCount));
557 prevDimCount = dimCount;
560 int totalCellCount = 1;
561 for (int dd = 0; dd < DIM; ++dd)
563 int cellCount;
564 if (bSingleCell[dd])
566 cellCount = 1;
568 else
570 cellCount = std::max(1, static_cast<int>(box[dd][dd] / targetsize));
571 // TODO: If the cell count is one or two, it would be better to
572 // just fall back to bSingleCell[dd] = true, and leave the rest to
573 // the efficiency check later.
574 if (bGridPBC_[dd] && cellCount < 3)
576 return false;
579 totalCellCount *= cellCount;
580 ncelldim_[dd] = cellCount;
582 if (totalCellCount <= 3)
584 return false;
586 // Never decrease the size of the cell vector to avoid reallocating
587 // memory for the nested vectors. The actual size of the vector is not
588 // used outside this function.
589 if (cells_.size() < static_cast<size_t>(totalCellCount))
591 cells_.resize(totalCellCount);
593 for (int ci = 0; ci < totalCellCount; ++ci)
595 cells_[ci].clear();
597 return true;
600 bool AnalysisNeighborhoodSearchImpl::initGrid(
601 const t_pbc &pbc, int posCount, const rvec x[], bool bForce)
603 if (posCount == 0)
605 return false;
608 switch (pbc.ePBC)
610 case epbcNONE:
611 bGridPBC_[XX] = false;
612 bGridPBC_[YY] = false;
613 bGridPBC_[ZZ] = false;
614 break;
615 case epbcXY:
616 bGridPBC_[XX] = true;
617 bGridPBC_[YY] = true;
618 bGridPBC_[ZZ] = false;
619 break;
620 case epbcXYZ:
621 bGridPBC_[XX] = true;
622 bGridPBC_[YY] = true;
623 bGridPBC_[ZZ] = true;
624 break;
625 default:
626 // Grid searching not supported for now with screw.
627 return false;
630 bool bSingleCell[DIM] = {false, false, bXY_};
631 matrix box;
632 copy_mat(pbc.box, box);
633 // TODO: In principle, we could use the bounding box for periodic
634 // dimensions as well if the bounding box is sufficiently far from the box
635 // edges.
636 rvec origin, boundingBoxSize;
637 computeBoundingBox(posCount, x, origin, boundingBoxSize);
638 clear_rvec(gridOrigin_);
639 for (int dd = 0; dd < DIM; ++dd)
641 if (!bGridPBC_[dd] && !bSingleCell[dd])
643 gridOrigin_[dd] = origin[dd];
644 clear_rvec(box[dd]);
645 box[dd][dd] = boundingBoxSize[dd];
647 // TODO: In case the zero vector comes from the bounding box, this does
648 // not lead to a very efficient grid search, but that should be rare.
649 if (box[dd][dd] <= 0.0)
651 GMX_ASSERT(!bGridPBC_[dd], "Periodic box vector is zero");
652 bSingleCell[dd] = true;
653 clear_rvec(box[dd]);
654 box[dd][dd] = 1.0;
658 if (!initGridCells(box, bSingleCell, posCount))
660 return false;
663 bTric_ = TRICLINIC(pbc.box);
664 for (int dd = 0; dd < DIM; ++dd)
666 cellSize_[dd] = box[dd][dd] / ncelldim_[dd];
667 if (bSingleCell[dd])
669 invCellSize_[dd] = 0.0;
671 else
673 invCellSize_[dd] = 1.0 / cellSize_[dd];
676 if (bTric_)
678 cellShiftZY_ = box[ZZ][YY] * invCellSize_[YY];
679 cellShiftZX_ = box[ZZ][XX] * invCellSize_[XX];
680 cellShiftYX_ = box[YY][XX] * invCellSize_[XX];
682 return checkGridSearchEfficiency(bForce);
685 void AnalysisNeighborhoodSearchImpl::mapPointToGridCell(const rvec x,
686 rvec cell,
687 rvec xout) const
689 rvec xtmp;
690 rvec_sub(x, gridOrigin_, xtmp);
691 // The reverse order is necessary for triclinic cells: shifting in Z may
692 // modify also X and Y, and shifting in Y may modify X, so the mapping to
693 // a rectangular grid needs to be done in this order.
694 for (int dd = DIM - 1; dd >= 0; --dd)
696 real cellIndex = xtmp[dd] * invCellSize_[dd];
697 if (bGridPBC_[dd])
699 const real cellCount = ncelldim_[dd];
700 while (cellIndex < 0)
702 cellIndex += cellCount;
703 rvec_inc(xtmp, pbc_.box[dd]);
705 while (cellIndex >= cellCount)
707 cellIndex -= cellCount;
708 rvec_dec(xtmp, pbc_.box[dd]);
711 cell[dd] = cellIndex;
713 copy_rvec(xtmp, xout);
716 int AnalysisNeighborhoodSearchImpl::getGridCellIndex(const ivec cell) const
718 GMX_ASSERT(cell[XX] >= 0 && cell[XX] < ncelldim_[XX],
719 "Grid cell X index out of range");
720 GMX_ASSERT(cell[YY] >= 0 && cell[YY] < ncelldim_[YY],
721 "Grid cell Y index out of range");
722 GMX_ASSERT(cell[ZZ] >= 0 && cell[ZZ] < ncelldim_[ZZ],
723 "Grid cell Z index out of range");
724 return cell[XX] + cell[YY] * ncelldim_[XX]
725 + cell[ZZ] * ncelldim_[XX] * ncelldim_[YY];
728 void AnalysisNeighborhoodSearchImpl::addToGridCell(const rvec cell, int i)
730 ivec icell;
731 for (int dd = 0; dd < DIM; ++dd)
733 int cellIndex = static_cast<int>(floor(cell[dd]));
734 if (!bGridPBC_[dd])
736 const int cellCount = ncelldim_[dd];
737 if (cellIndex < 0)
739 cellIndex = 0;
741 else if (cellIndex >= cellCount)
743 cellIndex = cellCount - 1;
746 icell[dd] = cellIndex;
748 const int ci = getGridCellIndex(icell);
749 cells_[ci].push_back(i);
752 void AnalysisNeighborhoodSearchImpl::initCellRange(
753 const rvec centerCell, ivec currCell, ivec upperBound, int dim) const
755 // TODO: Prune off cells that are completely outside the cutoff.
756 const real range = cutoff_ * invCellSize_[dim];
757 real startOffset = centerCell[dim] - range;
758 real endOffset = centerCell[dim] + range;
759 if (bTric_)
761 switch (dim)
763 case ZZ:
764 break;
765 case YY:
766 if (currCell[ZZ] < 0)
768 startOffset += cellShiftZY_;
769 endOffset += cellShiftZY_;
771 else if (currCell[ZZ] >= ncelldim_[ZZ])
773 startOffset -= cellShiftZY_;
774 endOffset -= cellShiftZY_;
776 break;
777 case XX:
778 if (currCell[ZZ] < 0)
780 startOffset += cellShiftZX_;
781 endOffset += cellShiftZX_;
783 else if (currCell[ZZ] >= ncelldim_[ZZ])
785 startOffset -= cellShiftZX_;
786 endOffset -= cellShiftZX_;
788 if (currCell[YY] < 0)
790 startOffset += cellShiftYX_;
791 endOffset += cellShiftYX_;
793 else if (currCell[YY] >= ncelldim_[YY])
795 startOffset -= cellShiftYX_;
796 endOffset -= cellShiftYX_;
798 break;
801 // For non-periodic dimensions, clamp to the actual grid edges.
802 if (!bGridPBC_[dim])
804 // If endOffset < 0 or startOffset > N, these may cause the whole
805 // test position/grid plane/grid row to be skipped.
806 if (startOffset < 0)
808 startOffset = 0;
810 const int cellCount = ncelldim_[dim];
811 if (endOffset > cellCount - 1)
813 endOffset = cellCount - 1;
816 currCell[dim] = static_cast<int>(floor(startOffset));
817 upperBound[dim] = static_cast<int>(floor(endOffset));
820 bool AnalysisNeighborhoodSearchImpl::nextCell(
821 const rvec centerCell, ivec cell, ivec upperBound) const
823 int dim = 0;
824 while (dim < DIM)
826 next:
827 ++cell[dim];
828 if (cell[dim] > upperBound[dim])
830 ++dim;
831 continue;
833 for (int d = dim - 1; d >= 0; --d)
835 initCellRange(centerCell, cell, upperBound, d);
836 if (cell[d] > upperBound[d])
838 dim = d + 1;
839 goto next;
842 return true;
844 return false;
847 int AnalysisNeighborhoodSearchImpl::shiftCell(const ivec cell, rvec shift) const
849 ivec shiftedCell;
850 copy_ivec(cell, shiftedCell);
852 clear_rvec(shift);
853 for (int d = 0; d < DIM; ++d)
855 const int cellCount = ncelldim_[d];
856 if (bGridPBC_[d])
858 // A single shift may not be sufficient if the cell must be shifted
859 // in more than one dimension, although for each individual
860 // dimension it would be.
861 while (shiftedCell[d] < 0)
863 shiftedCell[d] += cellCount;
864 rvec_inc(shift, pbc_.box[d]);
866 while (shiftedCell[d] >= cellCount)
868 shiftedCell[d] -= cellCount;
869 rvec_dec(shift, pbc_.box[d]);
874 return getGridCellIndex(shiftedCell);
877 void AnalysisNeighborhoodSearchImpl::init(
878 AnalysisNeighborhood::SearchMode mode,
879 bool bXY,
880 const t_blocka *excls,
881 const t_pbc *pbc,
882 const AnalysisNeighborhoodPositions &positions)
884 GMX_RELEASE_ASSERT(positions.index_ == -1,
885 "Individual indexed positions not supported as reference");
886 bXY_ = bXY;
887 if (bXY_ && pbc != nullptr && pbc->ePBC != epbcNONE)
889 if (pbc->ePBC != epbcXY && pbc->ePBC != epbcXYZ)
891 std::string message =
892 formatString("Computations in the XY plane are not supported with PBC type '%s'",
893 epbc_names[pbc->ePBC]);
894 GMX_THROW(NotImplementedError(message));
896 if (pbc->ePBC == epbcXYZ &&
897 (std::fabs(pbc->box[ZZ][XX]) > GMX_REAL_EPS*pbc->box[ZZ][ZZ] ||
898 std::fabs(pbc->box[ZZ][YY]) > GMX_REAL_EPS*pbc->box[ZZ][ZZ]))
900 GMX_THROW(NotImplementedError("Computations in the XY plane are not supported when the last box vector is not parallel to the Z axis"));
902 // Use a single grid cell in Z direction.
903 matrix box;
904 copy_mat(pbc->box, box);
905 clear_rvec(box[ZZ]);
906 set_pbc(&pbc_, epbcXY, box);
908 else if (pbc != nullptr)
910 pbc_ = *pbc;
912 else
914 pbc_.ePBC = epbcNONE;
915 clear_mat(pbc_.box);
917 nref_ = positions.count_;
918 if (mode == AnalysisNeighborhood::eSearchMode_Simple)
920 bGrid_ = false;
922 else if (bTryGrid_)
924 bGrid_ = initGrid(pbc_, positions.count_, positions.x_,
925 mode == AnalysisNeighborhood::eSearchMode_Grid);
927 refIndices_ = positions.indices_;
928 if (bGrid_)
930 xrefAlloc_.resize(nref_);
931 xref_ = as_rvec_array(xrefAlloc_.data());
933 for (int i = 0; i < nref_; ++i)
935 const int ii = (refIndices_ != nullptr) ? refIndices_[i] : i;
936 rvec refcell;
937 mapPointToGridCell(positions.x_[ii], refcell, xrefAlloc_[i]);
938 addToGridCell(refcell, i);
941 else if (refIndices_ != nullptr)
943 xrefAlloc_.resize(nref_);
944 xref_ = as_rvec_array(xrefAlloc_.data());
945 for (int i = 0; i < nref_; ++i)
947 copy_rvec(positions.x_[refIndices_[i]], xrefAlloc_[i]);
950 else
952 xref_ = positions.x_;
954 excls_ = excls;
955 refExclusionIds_ = nullptr;
956 if (excls != nullptr)
958 // TODO: Check that the IDs are ascending, or remove the limitation.
959 refExclusionIds_ = positions.exclusionIds_;
960 GMX_RELEASE_ASSERT(refExclusionIds_ != nullptr,
961 "Exclusion IDs must be set for reference positions "
962 "when exclusions are enabled");
966 /********************************************************************
967 * AnalysisNeighborhoodPairSearchImpl
970 void AnalysisNeighborhoodPairSearchImpl::reset(int testIndex)
972 testIndex_ = testIndex;
973 if (testIndex_ >= 0 && testIndex_ < testPosCount_)
975 const int index =
976 (testIndices_ != nullptr ? testIndices_[testIndex] : testIndex);
977 if (search_.bGrid_)
979 search_.mapPointToGridCell(testPositions_[index], testcell_, xtest_);
980 search_.initCellRange(testcell_, currCell_, cellBound_, ZZ);
981 search_.initCellRange(testcell_, currCell_, cellBound_, YY);
982 search_.initCellRange(testcell_, currCell_, cellBound_, XX);
984 else
986 copy_rvec(testPositions_[index], xtest_);
988 if (search_.excls_ != nullptr)
990 const int exclIndex = testExclusionIds_[index];
991 if (exclIndex < search_.excls_->nr)
993 const int startIndex = search_.excls_->index[exclIndex];
994 nexcl_ = search_.excls_->index[exclIndex + 1] - startIndex;
995 excl_ = &search_.excls_->a[startIndex];
997 else
999 nexcl_ = 0;
1000 excl_ = nullptr;
1004 previ_ = -1;
1005 prevr2_ = 0.0;
1006 clear_rvec(prevdx_);
1007 exclind_ = 0;
1008 prevcai_ = -1;
1011 void AnalysisNeighborhoodPairSearchImpl::nextTestPosition()
1013 if (testIndex_ < testPosCount_)
1015 ++testIndex_;
1016 reset(testIndex_);
1020 bool AnalysisNeighborhoodPairSearchImpl::isExcluded(int j)
1022 if (exclind_ < nexcl_)
1024 const int index =
1025 (search_.refIndices_ != nullptr ? search_.refIndices_[j] : j);
1026 const int refId = search_.refExclusionIds_[index];
1027 while (exclind_ < nexcl_ && excl_[exclind_] < refId)
1029 ++exclind_;
1031 if (exclind_ < nexcl_ && refId == excl_[exclind_])
1033 ++exclind_;
1034 return true;
1037 return false;
1040 void AnalysisNeighborhoodPairSearchImpl::startSearch(
1041 const AnalysisNeighborhoodPositions &positions)
1043 testPosCount_ = positions.count_;
1044 testPositions_ = positions.x_;
1045 testExclusionIds_ = positions.exclusionIds_;
1046 testIndices_ = positions.indices_;
1047 GMX_RELEASE_ASSERT(search_.excls_ == nullptr || testExclusionIds_ != nullptr,
1048 "Exclusion IDs must be set when exclusions are enabled");
1049 if (positions.index_ < 0)
1051 reset(0);
1053 else
1055 // Somewhat of a hack: setup the array such that only the last position
1056 // will be used.
1057 testPosCount_ = positions.index_ + 1;
1058 reset(positions.index_);
1062 template <class Action>
1063 bool AnalysisNeighborhoodPairSearchImpl::searchNext(Action action)
1065 while (testIndex_ < testPosCount_)
1067 if (search_.bGrid_)
1069 int cai = prevcai_ + 1;
1073 rvec shift;
1074 const int ci = search_.shiftCell(currCell_, shift);
1075 const int cellSize = static_cast<int>(search_.cells_[ci].size());
1076 for (; cai < cellSize; ++cai)
1078 const int i = search_.cells_[ci][cai];
1079 if (isExcluded(i))
1081 continue;
1083 rvec dx;
1084 rvec_sub(search_.xref_[i], xtest_, dx);
1085 rvec_sub(dx, shift, dx);
1086 const real r2
1087 = search_.bXY_
1088 ? dx[XX]*dx[XX] + dx[YY]*dx[YY]
1089 : norm2(dx);
1090 if (r2 <= search_.cutoff2_)
1092 if (action(i, r2, dx))
1094 prevcai_ = cai;
1095 previ_ = i;
1096 prevr2_ = r2;
1097 copy_rvec(dx, prevdx_);
1098 return true;
1102 exclind_ = 0;
1103 cai = 0;
1105 while (search_.nextCell(testcell_, currCell_, cellBound_));
1107 else
1109 for (int i = previ_ + 1; i < search_.nref_; ++i)
1111 if (isExcluded(i))
1113 continue;
1115 rvec dx;
1116 if (search_.pbc_.ePBC != epbcNONE)
1118 pbc_dx(&search_.pbc_, search_.xref_[i], xtest_, dx);
1120 else
1122 rvec_sub(search_.xref_[i], xtest_, dx);
1124 const real r2
1125 = search_.bXY_
1126 ? dx[XX]*dx[XX] + dx[YY]*dx[YY]
1127 : norm2(dx);
1128 if (r2 <= search_.cutoff2_)
1130 if (action(i, r2, dx))
1132 previ_ = i;
1133 prevr2_ = r2;
1134 copy_rvec(dx, prevdx_);
1135 return true;
1140 nextTestPosition();
1142 return false;
1145 void AnalysisNeighborhoodPairSearchImpl::initFoundPair(
1146 AnalysisNeighborhoodPair *pair) const
1148 if (previ_ < 0)
1150 *pair = AnalysisNeighborhoodPair();
1152 else
1154 *pair = AnalysisNeighborhoodPair(previ_, testIndex_, prevr2_, prevdx_);
1158 } // namespace internal
1160 namespace
1163 /*! \brief
1164 * Search action to find the next neighbor.
1166 * Used as the action for AnalysisNeighborhoodPairSearchImpl::searchNext() to
1167 * find the next neighbor.
1169 * Simply breaks the loop on the first found neighbor.
1171 bool withinAction(int /*i*/, real /*r2*/, const rvec /* dx */)
1173 return true;
1176 /*! \brief
1177 * Search action find the minimum distance.
1179 * Used as the action for AnalysisNeighborhoodPairSearchImpl::searchNext() to
1180 * find the nearest neighbor.
1182 * With this action, AnalysisNeighborhoodPairSearchImpl::searchNext() always
1183 * returns false, and the output is put into the variables passed by pointer
1184 * into the constructor. If no neighbors are found, the output variables are
1185 * not modified, i.e., the caller must initialize them.
1187 class MindistAction
1189 public:
1190 /*! \brief
1191 * Initializes the action with given output locations.
1193 * \param[out] closestPoint Index of the closest reference location.
1194 * \param[out] minDist2 Minimum distance squared.
1195 * \param[out] dx Shortest distance vector.
1197 * The constructor call does not modify the pointed values, but only
1198 * stores the pointers for later use.
1199 * See the class description for additional semantics.
1201 MindistAction(int *closestPoint, real *minDist2, rvec *dx)
1202 : closestPoint_(*closestPoint), minDist2_(*minDist2), dx_(*dx)
1205 //! Copies the action.
1206 MindistAction(const MindistAction &) = default;
1208 //! Processes a neighbor to find the nearest point.
1209 bool operator()(int i, real r2, const rvec dx)
1211 if (r2 < minDist2_)
1213 closestPoint_ = i;
1214 minDist2_ = r2;
1215 copy_rvec(dx, dx_);
1217 return false;
1220 private:
1221 int &closestPoint_;
1222 real &minDist2_;
1223 rvec &dx_;
1225 GMX_DISALLOW_ASSIGN(MindistAction);
1228 } // namespace
1230 /********************************************************************
1231 * AnalysisNeighborhood::Impl
1234 class AnalysisNeighborhood::Impl
1236 public:
1237 typedef AnalysisNeighborhoodSearch::ImplPointer SearchImplPointer;
1238 typedef std::vector<SearchImplPointer> SearchList;
1240 Impl()
1241 : cutoff_(0), excls_(nullptr), mode_(eSearchMode_Automatic), bXY_(false)
1244 ~Impl()
1246 SearchList::const_iterator i;
1247 for (i = searchList_.begin(); i != searchList_.end(); ++i)
1249 GMX_RELEASE_ASSERT(i->unique(),
1250 "Dangling AnalysisNeighborhoodSearch reference");
1254 SearchImplPointer getSearch();
1256 Mutex createSearchMutex_;
1257 SearchList searchList_;
1258 real cutoff_;
1259 const t_blocka *excls_;
1260 SearchMode mode_;
1261 bool bXY_;
1264 AnalysisNeighborhood::Impl::SearchImplPointer
1265 AnalysisNeighborhood::Impl::getSearch()
1267 lock_guard<Mutex> lock(createSearchMutex_);
1268 // TODO: Consider whether this needs to/can be faster, e.g., by keeping a
1269 // separate pool of unused search objects.
1270 SearchList::const_iterator i;
1271 for (i = searchList_.begin(); i != searchList_.end(); ++i)
1273 if (i->unique())
1275 return *i;
1278 SearchImplPointer search(new internal::AnalysisNeighborhoodSearchImpl(cutoff_));
1279 searchList_.push_back(search);
1280 return search;
1283 /********************************************************************
1284 * AnalysisNeighborhood
1287 AnalysisNeighborhood::AnalysisNeighborhood()
1288 : impl_(new Impl)
1292 AnalysisNeighborhood::~AnalysisNeighborhood()
1296 void AnalysisNeighborhood::setCutoff(real cutoff)
1298 GMX_RELEASE_ASSERT(impl_->searchList_.empty(),
1299 "Changing the cutoff after initSearch() not currently supported");
1300 impl_->cutoff_ = cutoff;
1303 void AnalysisNeighborhood::setXYMode(bool bXY)
1305 impl_->bXY_ = bXY;
1308 void AnalysisNeighborhood::setTopologyExclusions(const t_blocka *excls)
1310 GMX_RELEASE_ASSERT(impl_->searchList_.empty(),
1311 "Changing the exclusions after initSearch() not currently supported");
1312 impl_->excls_ = excls;
1315 void AnalysisNeighborhood::setMode(SearchMode mode)
1317 impl_->mode_ = mode;
1320 AnalysisNeighborhood::SearchMode AnalysisNeighborhood::mode() const
1322 return impl_->mode_;
1325 AnalysisNeighborhoodSearch
1326 AnalysisNeighborhood::initSearch(const t_pbc *pbc,
1327 const AnalysisNeighborhoodPositions &positions)
1329 Impl::SearchImplPointer search(impl_->getSearch());
1330 search->init(mode(), impl_->bXY_, impl_->excls_,
1331 pbc, positions);
1332 return AnalysisNeighborhoodSearch(search);
1335 /********************************************************************
1336 * AnalysisNeighborhoodSearch
1339 AnalysisNeighborhoodSearch::AnalysisNeighborhoodSearch()
1343 AnalysisNeighborhoodSearch::AnalysisNeighborhoodSearch(const ImplPointer &impl)
1344 : impl_(impl)
1348 void AnalysisNeighborhoodSearch::reset()
1350 impl_.reset();
1353 AnalysisNeighborhood::SearchMode AnalysisNeighborhoodSearch::mode() const
1355 GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
1356 return (impl_->usesGridSearch()
1357 ? AnalysisNeighborhood::eSearchMode_Grid
1358 : AnalysisNeighborhood::eSearchMode_Simple);
1361 bool AnalysisNeighborhoodSearch::isWithin(
1362 const AnalysisNeighborhoodPositions &positions) const
1364 GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
1365 internal::AnalysisNeighborhoodPairSearchImpl pairSearch(*impl_);
1366 pairSearch.startSearch(positions);
1367 return pairSearch.searchNext(&withinAction);
1370 real AnalysisNeighborhoodSearch::minimumDistance(
1371 const AnalysisNeighborhoodPositions &positions) const
1373 GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
1374 internal::AnalysisNeighborhoodPairSearchImpl pairSearch(*impl_);
1375 pairSearch.startSearch(positions);
1376 real minDist2 = impl_->cutoffSquared();
1377 int closestPoint = -1;
1378 rvec dx = {0.0, 0.0, 0.0};
1379 MindistAction action(&closestPoint, &minDist2, &dx);
1380 (void)pairSearch.searchNext(action);
1381 return std::sqrt(minDist2);
1384 AnalysisNeighborhoodPair
1385 AnalysisNeighborhoodSearch::nearestPoint(
1386 const AnalysisNeighborhoodPositions &positions) const
1388 GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
1389 internal::AnalysisNeighborhoodPairSearchImpl pairSearch(*impl_);
1390 pairSearch.startSearch(positions);
1391 real minDist2 = impl_->cutoffSquared();
1392 int closestPoint = -1;
1393 rvec dx = {0.0, 0.0, 0.0};
1394 MindistAction action(&closestPoint, &minDist2, &dx);
1395 (void)pairSearch.searchNext(action);
1396 return AnalysisNeighborhoodPair(closestPoint, 0, minDist2, dx);
1399 AnalysisNeighborhoodPairSearch
1400 AnalysisNeighborhoodSearch::startPairSearch(
1401 const AnalysisNeighborhoodPositions &positions) const
1403 GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
1404 Impl::PairSearchImplPointer pairSearch(impl_->getPairSearch());
1405 pairSearch->startSearch(positions);
1406 return AnalysisNeighborhoodPairSearch(pairSearch);
1409 /********************************************************************
1410 * AnalysisNeighborhoodPairSearch
1413 AnalysisNeighborhoodPairSearch::AnalysisNeighborhoodPairSearch(
1414 const ImplPointer &impl)
1415 : impl_(impl)
1419 bool AnalysisNeighborhoodPairSearch::findNextPair(AnalysisNeighborhoodPair *pair)
1421 bool bFound = impl_->searchNext(&withinAction);
1422 impl_->initFoundPair(pair);
1423 return bFound;
1426 void AnalysisNeighborhoodPairSearch::skipRemainingPairsForTestPosition()
1428 impl_->nextTestPosition();
1431 } // namespace gmx