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.
37 * Implements neighborhood searching for analysis (from nbsearch.h).
39 * High-level overview of the algorithm is at \ref page_analysisnbsearch.
42 * The grid implementation could still be optimized in several different ways:
43 * - A better heuristic for selecting the grid size or falling back to a
44 * simple all-pairs search.
45 * - A multi-level grid implementation could be used to be able to use small
46 * grids for short cutoffs with very inhomogeneous particle distributions
47 * without a memory cost.
49 * \author Teemu Murtola <teemu.murtola@gmail.com>
50 * \ingroup module_selection
62 #include "gromacs/math/functions.h"
63 #include "gromacs/math/vec.h"
64 #include "gromacs/pbcutil/pbc.h"
65 #include "gromacs/selection/position.h"
66 #include "gromacs/topology/block.h"
67 #include "gromacs/utility/arrayref.h"
68 #include "gromacs/utility/exceptions.h"
69 #include "gromacs/utility/gmxassert.h"
70 #include "gromacs/utility/mutex.h"
71 #include "gromacs/utility/stringutil.h"
80 * Computes the bounding box for a set of positions.
82 * \param[in] posCount Number of positions in \p x.
83 * \param[in] x Positions to compute the bounding box for.
84 * \param[out] origin Origin of the bounding box.
85 * \param[out] size Size of the bounding box.
87 void computeBoundingBox(int posCount
, const rvec x
[], rvec origin
, rvec size
)
90 copy_rvec(x
[0], origin
);
91 copy_rvec(x
[0], maxBound
);
92 for (int i
= 1; i
< posCount
; ++i
)
94 for (int d
= 0; d
< DIM
; ++d
)
96 if (origin
[d
] > x
[i
][d
])
100 if (maxBound
[d
] < x
[i
][d
])
102 maxBound
[d
] = x
[i
][d
];
106 rvec_sub(maxBound
, origin
, size
);
114 /********************************************************************
115 * Implementation class declarations
118 class AnalysisNeighborhoodSearchImpl
121 typedef AnalysisNeighborhoodPairSearch::ImplPointer
122 PairSearchImplPointer
;
123 typedef std::vector
<PairSearchImplPointer
> PairSearchList
;
124 typedef std::vector
<std::vector
<int> > CellList
;
126 explicit AnalysisNeighborhoodSearchImpl(real cutoff
);
127 ~AnalysisNeighborhoodSearchImpl();
130 * Initializes the search with a given box and reference positions.
132 * \param[in] mode Search mode to use.
133 * \param[in] bXY Whether to use 2D searching.
134 * \param[in] excls Exclusions.
135 * \param[in] pbc PBC information.
136 * \param[in] positions Set of reference positions.
138 void init(AnalysisNeighborhood::SearchMode mode
,
140 const t_blocka
*excls
,
142 const AnalysisNeighborhoodPositions
&positions
);
143 PairSearchImplPointer
getPairSearch();
145 real
cutoffSquared() const { return cutoff2_
; }
146 bool usesGridSearch() const { return bGrid_
; }
150 * Determines a suitable grid size and sets up the cells.
152 * \param[in] box Box vectors (should not have zero vectors).
153 * \param[in] bSingleCell If `true`, the corresponding dimension will
154 * be forced to use a single cell.
155 * \param[in] posCount Number of positions that will be put on the
157 * \returns `false` if grid search is not suitable.
159 bool initGridCells(const matrix box
, bool bSingleCell
[DIM
],
162 * Sets ua a search grid for a given box.
164 * \param[in] pbc Information about the box.
165 * \param[in] posCount Number of positions in \p x.
166 * \param[in] x Reference positions that will be put on the grid.
167 * \param[in] bForce If `true`, grid searching will be used if at all
168 * possible, even if a simple search might give better performance.
169 * \returns `false` if grid search is not suitable.
171 bool initGrid(const t_pbc
&pbc
, int posCount
, const rvec x
[], bool bForce
);
173 * Maps a point into a grid cell.
175 * \param[in] x Point to map.
176 * \param[out] cell Fractional cell coordinates of \p x on the grid.
177 * \param[out] xout Coordinates to use.
179 * \p xout will be within the rectangular unit cell in dimensions where
180 * the grid is periodic. For other dimensions, both \p xout and
181 * \p cell can be outside the grid/unit cell.
183 void mapPointToGridCell(const rvec x
, rvec cell
, rvec xout
) const;
185 * Calculates linear index of a grid cell.
187 * \param[in] cell Cell indices (must be within the grid).
188 * \returns Linear index of \p cell.
190 int getGridCellIndex(const ivec cell
) const;
192 * Calculates linear index of a grid cell from fractional coordinates.
194 * \param[in] cell Cell indices (must be within the grid).
195 * \returns Linear index of \p cell.
197 int getGridCellIndex(const rvec cell
) const;
199 * Adds an index into a grid cell.
201 * \param[in] cell Fractional cell coordinates into which \p i should
203 * \param[in] i Index to add.
205 * \p cell should satisfy the conditions that \p mapPointToGridCell()
208 void addToGridCell(const rvec cell
, int i
);
210 * Initializes a cell pair loop for a dimension.
212 * \param[in] centerCell Fractional cell coordiates of the particle
213 * for which pairs are being searched.
214 * \param[in,out] cell Current/initial cell to loop over.
215 * \param[in,out] upperBound Last cell to loop over.
216 * \param[in] dim Dimension to initialize in this call.
218 * Initializes `cell[dim]` and `upperBound[dim]` for looping over
219 * neighbors of a particle at position given by \p centerCell.
220 * If 'dim != ZZ`, `cell[d]` (`d > dim`) set the plane/row of cells
221 * for which the loop is initialized. The loop should then go from
222 * `cell[dim]` until `upperBound[dim]`, inclusive.
223 * `cell[d]` with `d < dim` or `upperBound[d]` with `d != dim` are not
224 * modified by this function.
226 * `cell` and `upperBound` may be outside the grid for periodic
227 * dimensions and need to be shifted separately: to simplify the
228 * looping, the range is always (roughly) symmetric around the value in
231 void initCellRange(const rvec centerCell
, ivec cell
,
232 ivec upperBound
, int dim
) const;
234 * Computes the extent of the cutoff sphere on a particular cell edge.
236 * \param[in] centerCell Fractional cell coordiates of the particle
237 * for which pairs are being searched.
238 * \param[in] cell Current cell (for dimensions `>dim`).
239 * \param[in] dim Dimension to compute in this call.
240 * \returns Fractional extent of the cutoff sphere when looping
241 * over cells in dimension `dim`, for `cell[d]` (`d > dim`).
243 * Input parameters are as for initCellRange(), except that if `cell`
244 * is over a periodic boundary from `centerCell`, triclinic shifts
245 * should have been applied to `centerCell` X/Y components.
247 real
computeCutoffExtent(const RVec centerCell
, const ivec cell
, int dim
) const;
249 * Advances cell pair loop to the next cell.
251 * \param[in] centerCell Fractional cell coordiates of the particle
252 * for which pairs are being searched.
253 * \param[in,out] cell Current (in)/next (out) cell in the loop.
254 * \param[in,out] upperBound Last cell in the loop for each dimension.
256 bool nextCell(const rvec centerCell
, ivec cell
, ivec upperBound
) const;
258 * Calculates the index and shift of a grid cell during looping.
260 * \param[in] cell Unshifted cell index.
261 * \param[out] shift Shift to apply to get the periodic distance
262 * for distances between the cells.
263 * \returns Grid cell index corresponding to `cell`.
265 int shiftCell(const ivec cell
, rvec shift
) const;
267 //! Whether to try grid searching.
271 //! The cutoff squared.
273 //! Whether to do searching in XY plane only.
276 //! Number of reference points for the current frame.
278 //! Reference point positions.
280 //! Reference position exclusion IDs.
281 const int *refExclusionIds_
;
282 //! Reference position indices (NULL if no indices).
283 const int *refIndices_
;
285 const t_blocka
*excls_
;
289 //! Whether grid searching is actually used for the current positions.
291 //! false if the box is rectangular.
293 //! Whether the grid is periodic in a dimension.
295 //! Array for storing in-unit-cell reference positions.
296 std::vector
<RVec
> xrefAlloc_
;
297 //! Origin of the grid (zero for periodic dimensions).
299 //! Size of a single grid cell.
301 //! Inverse of \p cellSize_. Zero for dimensions where grid is not used.
304 * Shift in cell coordinates (for triclinic boxes) in X when crossing
305 * the Z periodic boundary.
309 * Shift in cell coordinates (for triclinic boxes) in Y when crossing
310 * the Z periodic boundary.
314 * Shift in cell coordinates (for triclinic boxes) in X when crossing
315 * the Y periodic boundary.
318 //! Number of cells along each dimension.
320 //! Data structure to hold the grid cell contents.
323 Mutex createPairSearchMutex_
;
324 PairSearchList pairSearchList_
;
326 friend class AnalysisNeighborhoodPairSearchImpl
;
328 GMX_DISALLOW_COPY_AND_ASSIGN(AnalysisNeighborhoodSearchImpl
);
331 class AnalysisNeighborhoodPairSearchImpl
334 explicit AnalysisNeighborhoodPairSearchImpl(const AnalysisNeighborhoodSearchImpl
&search
)
337 selfSearchMode_
= false;
339 testPositions_
= nullptr;
340 testExclusionIds_
= nullptr;
341 testIndices_
= nullptr;
345 clear_rvec(testcell_
);
346 clear_ivec(currCell_
);
347 clear_ivec(cellBound_
);
351 //! Initializes a search to find reference positions neighboring \p x.
352 void startSearch(const AnalysisNeighborhoodPositions
&positions
);
353 //! Initializes a search to find reference position pairs.
354 void startSelfSearch();
355 //! Searches for the next neighbor.
356 template <class Action
>
357 bool searchNext(Action action
);
358 //! Initializes a pair representing the pair found by searchNext().
359 void initFoundPair(AnalysisNeighborhoodPair
*pair
) const;
360 //! Advances to the next test position, skipping any remaining pairs.
361 void nextTestPosition();
364 //! Clears the loop indices.
365 void reset(int testIndex
);
366 //! Checks whether a reference positiong should be excluded.
367 bool isExcluded(int j
);
369 //! Parent search object.
370 const AnalysisNeighborhoodSearchImpl
&search_
;
371 //! Whether we are searching for ref-ref pairs.
372 bool selfSearchMode_
;
373 //! Number of test positions.
375 //! Reference to the test positions.
376 const rvec
*testPositions_
;
377 //! Reference to the test exclusion indices.
378 const int *testExclusionIds_
;
379 //! Reference to the test position indices.
380 const int *testIndices_
;
381 //! Number of excluded reference positions for current test particle.
383 //! Exclusions for current test particle.
385 //! Index of the currently active test position in \p testPositions_.
387 //! Stores test position during a pair loop.
389 //! Stores the previous returned position during a pair loop.
391 //! Stores the pair distance corresponding to previ_.
393 //! Stores the shortest distance vector corresponding to previ_.
395 //! Stores the current exclusion index during loops.
397 //! Stores the fractional test particle cell location during loops.
399 //! Stores the cell index corresponding to testcell_.
401 //! Stores the current cell during pair loops.
403 //! Stores the current loop upper bounds for each dimension during pair loops.
405 //! Stores the index within the current cell during pair loops.
408 GMX_DISALLOW_COPY_AND_ASSIGN(AnalysisNeighborhoodPairSearchImpl
);
411 /********************************************************************
412 * AnalysisNeighborhoodSearchImpl
415 AnalysisNeighborhoodSearchImpl::AnalysisNeighborhoodSearchImpl(real cutoff
)
421 cutoff_
= cutoff2_
= GMX_REAL_MAX
;
426 cutoff2_
= gmx::square(cutoff_
);
431 refExclusionIds_
= nullptr;
432 refIndices_
= nullptr;
433 std::memset(&pbc_
, 0, sizeof(pbc_
));
437 bGridPBC_
[XX
] = true;
438 bGridPBC_
[YY
] = true;
439 bGridPBC_
[ZZ
] = true;
441 clear_rvec(gridOrigin_
);
442 clear_rvec(cellSize_
);
443 clear_rvec(invCellSize_
);
444 clear_ivec(ncelldim_
);
447 AnalysisNeighborhoodSearchImpl::~AnalysisNeighborhoodSearchImpl()
449 PairSearchList::const_iterator i
;
450 for (i
= pairSearchList_
.begin(); i
!= pairSearchList_
.end(); ++i
)
452 GMX_RELEASE_ASSERT(i
->unique(),
453 "Dangling AnalysisNeighborhoodPairSearch reference");
457 AnalysisNeighborhoodSearchImpl::PairSearchImplPointer
458 AnalysisNeighborhoodSearchImpl::getPairSearch()
460 lock_guard
<Mutex
> lock(createPairSearchMutex_
);
461 // TODO: Consider whether this needs to/can be faster, e.g., by keeping a
462 // separate pool of unused search objects.
463 PairSearchList::const_iterator i
;
464 for (i
= pairSearchList_
.begin(); i
!= pairSearchList_
.end(); ++i
)
471 PairSearchImplPointer
pairSearch(new AnalysisNeighborhoodPairSearchImpl(*this));
472 pairSearchList_
.push_back(pairSearch
);
476 bool AnalysisNeighborhoodSearchImpl::initGridCells(
477 const matrix box
, bool bSingleCell
[DIM
], int posCount
)
479 // Determine the size of cubes where there are on average 10 positions.
480 // The loop takes care of cases where some of the box edges are shorter
481 // than the the desired cube size; in such cases, a single grid cell is
482 // used in these dimensions, and the cube size is determined only from the
483 // larger box vectors. Such boxes should be rare, but the bounding box
484 // approach can result in very flat boxes with certain types of selections
485 // (e.g., for interfacial systems or for small number of atoms).
486 real targetsize
= 0.0;
487 int prevDimCount
= 4;
492 for (int dd
= 0; dd
< DIM
; ++dd
)
494 const real boxSize
= box
[dd
][dd
];
495 if (boxSize
< targetsize
)
497 bSingleCell
[dd
] = true;
500 // TODO: Consider if a fallback would be possible/better.
513 if (dimCount
== 0 || dimCount
== prevDimCount
)
517 targetsize
= pow(volume
* 10 / posCount
, static_cast<real
>(1./dimCount
));
518 prevDimCount
= dimCount
;
521 int totalCellCount
= 1;
522 for (int dd
= 0; dd
< DIM
; ++dd
)
531 cellCount
= std::max(1, static_cast<int>(box
[dd
][dd
] / targetsize
));
532 // TODO: If the cell count is one or two, it could be better to
533 // just fall back to bSingleCell[dd] = true.
534 if (bGridPBC_
[dd
] && cellCount
< 3)
539 totalCellCount
*= cellCount
;
540 ncelldim_
[dd
] = cellCount
;
542 if (totalCellCount
<= 3)
546 // Never decrease the size of the cell vector to avoid reallocating
547 // memory for the nested vectors. The actual size of the vector is not
548 // used outside this function.
549 if (cells_
.size() < static_cast<size_t>(totalCellCount
))
551 cells_
.resize(totalCellCount
);
553 for (int ci
= 0; ci
< totalCellCount
; ++ci
)
560 bool AnalysisNeighborhoodSearchImpl::initGrid(
561 const t_pbc
&pbc
, int posCount
, const rvec x
[], bool bForce
)
568 // TODO: Use this again (can be useful when tuning initGridCells()),
569 // or remove throughout.
570 GMX_UNUSED_VALUE(bForce
);
575 bGridPBC_
[XX
] = false;
576 bGridPBC_
[YY
] = false;
577 bGridPBC_
[ZZ
] = false;
580 bGridPBC_
[XX
] = true;
581 bGridPBC_
[YY
] = true;
582 bGridPBC_
[ZZ
] = false;
585 bGridPBC_
[XX
] = true;
586 bGridPBC_
[YY
] = true;
587 bGridPBC_
[ZZ
] = true;
590 // Grid searching not supported for now with screw.
594 bool bSingleCell
[DIM
] = {false, false, bXY_
};
596 copy_mat(pbc
.box
, box
);
597 // TODO: In principle, we could use the bounding box for periodic
598 // dimensions as well if the bounding box is sufficiently far from the box
600 rvec origin
, boundingBoxSize
;
601 computeBoundingBox(posCount
, x
, origin
, boundingBoxSize
);
602 clear_rvec(gridOrigin_
);
603 for (int dd
= 0; dd
< DIM
; ++dd
)
605 if (!bGridPBC_
[dd
] && !bSingleCell
[dd
])
607 gridOrigin_
[dd
] = origin
[dd
];
609 box
[dd
][dd
] = boundingBoxSize
[dd
];
611 // TODO: In case the zero vector comes from the bounding box, this does
612 // not lead to a very efficient grid search, but that should be rare.
613 if (box
[dd
][dd
] <= 0.0)
615 GMX_ASSERT(!bGridPBC_
[dd
], "Periodic box vector is zero");
616 bSingleCell
[dd
] = true;
622 if (!initGridCells(box
, bSingleCell
, posCount
))
627 bTric_
= TRICLINIC(pbc
.box
);
628 for (int dd
= 0; dd
< DIM
; ++dd
)
630 cellSize_
[dd
] = box
[dd
][dd
] / ncelldim_
[dd
];
633 invCellSize_
[dd
] = 0.0;
637 invCellSize_
[dd
] = 1.0 / cellSize_
[dd
];
638 // TODO: It could be better to avoid this when determining the cell
639 // size, but this can still remain here as a fallback to avoid
640 // incorrect results.
641 if (std::ceil(2*cutoff_
*invCellSize_
[dd
]) >= ncelldim_
[dd
])
643 // Cutoff is too close to half the box size for grid searching
644 // (it is not possible to find a single shift for every pair of
652 cellShiftZY_
= box
[ZZ
][YY
] * invCellSize_
[YY
];
653 cellShiftZX_
= box
[ZZ
][XX
] * invCellSize_
[XX
];
654 cellShiftYX_
= box
[YY
][XX
] * invCellSize_
[XX
];
659 void AnalysisNeighborhoodSearchImpl::mapPointToGridCell(const rvec x
,
664 rvec_sub(x
, gridOrigin_
, xtmp
);
665 // The reverse order is necessary for triclinic cells: shifting in Z may
666 // modify also X and Y, and shifting in Y may modify X, so the mapping to
667 // a rectangular grid needs to be done in this order.
668 for (int dd
= DIM
- 1; dd
>= 0; --dd
)
670 real cellIndex
= xtmp
[dd
] * invCellSize_
[dd
];
673 const real cellCount
= ncelldim_
[dd
];
674 while (cellIndex
< 0)
676 cellIndex
+= cellCount
;
677 rvec_inc(xtmp
, pbc_
.box
[dd
]);
679 while (cellIndex
>= cellCount
)
681 cellIndex
-= cellCount
;
682 rvec_dec(xtmp
, pbc_
.box
[dd
]);
685 cell
[dd
] = cellIndex
;
687 copy_rvec(xtmp
, xout
);
690 int AnalysisNeighborhoodSearchImpl::getGridCellIndex(const ivec cell
) const
692 GMX_ASSERT(cell
[XX
] >= 0 && cell
[XX
] < ncelldim_
[XX
],
693 "Grid cell X index out of range");
694 GMX_ASSERT(cell
[YY
] >= 0 && cell
[YY
] < ncelldim_
[YY
],
695 "Grid cell Y index out of range");
696 GMX_ASSERT(cell
[ZZ
] >= 0 && cell
[ZZ
] < ncelldim_
[ZZ
],
697 "Grid cell Z index out of range");
698 return cell
[XX
] + cell
[YY
] * ncelldim_
[XX
]
699 + cell
[ZZ
] * ncelldim_
[XX
] * ncelldim_
[YY
];
702 int AnalysisNeighborhoodSearchImpl::getGridCellIndex(const rvec cell
) const
705 for (int dd
= 0; dd
< DIM
; ++dd
)
707 int cellIndex
= static_cast<int>(floor(cell
[dd
]));
710 const int cellCount
= ncelldim_
[dd
];
715 else if (cellIndex
>= cellCount
)
717 cellIndex
= cellCount
- 1;
720 icell
[dd
] = cellIndex
;
722 return getGridCellIndex(icell
);
725 void AnalysisNeighborhoodSearchImpl::addToGridCell(const rvec cell
, int i
)
727 const int ci
= getGridCellIndex(cell
);
728 cells_
[ci
].push_back(i
);
731 void AnalysisNeighborhoodSearchImpl::initCellRange(
732 const rvec centerCell
, ivec currCell
, ivec upperBound
, int dim
) const
734 RVec
shiftedCenter(centerCell
);
735 // Shift the center to the cell coordinates of currCell, so that
736 // computeCutoffExtent() can assume simple rectangular grid.
741 if (currCell
[ZZ
] < 0)
743 shiftedCenter
[XX
] += cellShiftZX_
;
745 else if (currCell
[ZZ
] >= ncelldim_
[ZZ
])
747 shiftedCenter
[XX
] -= cellShiftZX_
;
749 if (currCell
[YY
] < 0)
751 shiftedCenter
[XX
] += cellShiftYX_
;
753 else if (currCell
[YY
] >= ncelldim_
[YY
])
755 shiftedCenter
[XX
] -= cellShiftYX_
;
758 if (dim
== XX
|| dim
== YY
)
760 if (currCell
[ZZ
] < 0)
762 shiftedCenter
[YY
] += cellShiftZY_
;
764 else if (currCell
[ZZ
] >= ncelldim_
[ZZ
])
766 shiftedCenter
[YY
] -= cellShiftZY_
;
770 const real range
= computeCutoffExtent(shiftedCenter
, currCell
, dim
) * invCellSize_
[dim
];
771 real startOffset
= shiftedCenter
[dim
] - range
;
772 real endOffset
= shiftedCenter
[dim
] + range
;
773 // For non-periodic dimensions, clamp to the actual grid edges.
776 // If endOffset < 0 or startOffset > N, these may cause the whole
777 // test position/grid plane/grid row to be skipped.
782 const int cellCount
= ncelldim_
[dim
];
783 if (endOffset
> cellCount
- 1)
785 endOffset
= cellCount
- 1;
788 currCell
[dim
] = static_cast<int>(floor(startOffset
));
789 upperBound
[dim
] = static_cast<int>(floor(endOffset
));
792 real
AnalysisNeighborhoodSearchImpl::computeCutoffExtent(
793 const RVec centerCell
, const ivec cell
, int dim
) const
801 for (int d
= dim
+ 1; d
< DIM
; ++d
)
803 real dimDist
= cell
[d
] - centerCell
[d
];
808 else if (dimDist
<= 0)
812 dist2
+= dimDist
*dimDist
*cellSize_
[d
]*cellSize_
[d
];
814 if (dist2
>= cutoff2_
)
818 return std::sqrt(cutoff2_
- dist2
);
821 bool AnalysisNeighborhoodSearchImpl::nextCell(
822 const rvec centerCell
, ivec cell
, ivec upperBound
) const
829 if (cell
[dim
] > upperBound
[dim
])
834 for (int d
= dim
- 1; d
>= 0; --d
)
836 initCellRange(centerCell
, cell
, upperBound
, d
);
837 if (cell
[d
] > upperBound
[d
])
848 int AnalysisNeighborhoodSearchImpl::shiftCell(const ivec cell
, rvec shift
) const
851 copy_ivec(cell
, shiftedCell
);
854 for (int d
= 0; d
< DIM
; ++d
)
856 const int cellCount
= ncelldim_
[d
];
859 // A single shift may not be sufficient if the cell must be shifted
860 // in more than one dimension, although for each individual
861 // dimension it would be.
862 while (shiftedCell
[d
] < 0)
864 shiftedCell
[d
] += cellCount
;
865 rvec_inc(shift
, pbc_
.box
[d
]);
867 while (shiftedCell
[d
] >= cellCount
)
869 shiftedCell
[d
] -= cellCount
;
870 rvec_dec(shift
, pbc_
.box
[d
]);
875 return getGridCellIndex(shiftedCell
);
878 void AnalysisNeighborhoodSearchImpl::init(
879 AnalysisNeighborhood::SearchMode mode
,
881 const t_blocka
*excls
,
883 const AnalysisNeighborhoodPositions
&positions
)
885 GMX_RELEASE_ASSERT(positions
.index_
== -1,
886 "Individual indexed positions not supported as reference");
888 if (bXY_
&& pbc
!= nullptr && pbc
->ePBC
!= epbcNONE
)
890 if (pbc
->ePBC
!= epbcXY
&& pbc
->ePBC
!= epbcXYZ
)
892 std::string message
=
893 formatString("Computations in the XY plane are not supported with PBC type '%s'",
894 epbc_names
[pbc
->ePBC
]);
895 GMX_THROW(NotImplementedError(message
));
897 if (pbc
->ePBC
== epbcXYZ
&&
898 (std::fabs(pbc
->box
[ZZ
][XX
]) > GMX_REAL_EPS
*pbc
->box
[ZZ
][ZZ
] ||
899 std::fabs(pbc
->box
[ZZ
][YY
]) > GMX_REAL_EPS
*pbc
->box
[ZZ
][ZZ
]))
901 GMX_THROW(NotImplementedError("Computations in the XY plane are not supported when the last box vector is not parallel to the Z axis"));
903 // Use a single grid cell in Z direction.
905 copy_mat(pbc
->box
, box
);
907 set_pbc(&pbc_
, epbcXY
, box
);
909 else if (pbc
!= nullptr)
915 pbc_
.ePBC
= epbcNONE
;
918 nref_
= positions
.count_
;
919 if (mode
== AnalysisNeighborhood::eSearchMode_Simple
)
925 bGrid_
= initGrid(pbc_
, positions
.count_
, positions
.x_
,
926 mode
== AnalysisNeighborhood::eSearchMode_Grid
);
928 refIndices_
= positions
.indices_
;
931 xrefAlloc_
.resize(nref_
);
932 xref_
= as_rvec_array(xrefAlloc_
.data());
934 for (int i
= 0; i
< nref_
; ++i
)
936 const int ii
= (refIndices_
!= nullptr) ? refIndices_
[i
] : i
;
938 mapPointToGridCell(positions
.x_
[ii
], refcell
, xrefAlloc_
[i
]);
939 addToGridCell(refcell
, i
);
942 else if (refIndices_
!= nullptr)
944 xrefAlloc_
.resize(nref_
);
945 xref_
= as_rvec_array(xrefAlloc_
.data());
946 for (int i
= 0; i
< nref_
; ++i
)
948 copy_rvec(positions
.x_
[refIndices_
[i
]], xrefAlloc_
[i
]);
953 xref_
= positions
.x_
;
956 refExclusionIds_
= nullptr;
957 if (excls
!= nullptr)
959 // TODO: Check that the IDs are ascending, or remove the limitation.
960 refExclusionIds_
= positions
.exclusionIds_
;
961 GMX_RELEASE_ASSERT(refExclusionIds_
!= nullptr,
962 "Exclusion IDs must be set for reference positions "
963 "when exclusions are enabled");
967 /********************************************************************
968 * AnalysisNeighborhoodPairSearchImpl
971 void AnalysisNeighborhoodPairSearchImpl::reset(int testIndex
)
973 testIndex_
= testIndex
;
980 if (testIndex_
>= 0 && testIndex_
< testPosCount_
)
983 (testIndices_
!= nullptr ? testIndices_
[testIndex
] : testIndex
);
986 search_
.mapPointToGridCell(testPositions_
[index
], testcell_
, xtest_
);
987 search_
.initCellRange(testcell_
, currCell_
, cellBound_
, ZZ
);
988 search_
.initCellRange(testcell_
, currCell_
, cellBound_
, YY
);
989 search_
.initCellRange(testcell_
, currCell_
, cellBound_
, XX
);
992 testCellIndex_
= search_
.getGridCellIndex(testcell_
);
997 copy_rvec(testPositions_
[index
], xtest_
);
1000 previ_
= testIndex_
;
1003 if (search_
.excls_
!= nullptr)
1005 const int exclIndex
= testExclusionIds_
[index
];
1006 if (exclIndex
< search_
.excls_
->nr
)
1008 const int startIndex
= search_
.excls_
->index
[exclIndex
];
1009 nexcl_
= search_
.excls_
->index
[exclIndex
+ 1] - startIndex
;
1010 excl_
= &search_
.excls_
->a
[startIndex
];
1021 void AnalysisNeighborhoodPairSearchImpl::nextTestPosition()
1023 if (testIndex_
< testPosCount_
)
1030 bool AnalysisNeighborhoodPairSearchImpl::isExcluded(int j
)
1032 if (exclind_
< nexcl_
)
1035 (search_
.refIndices_
!= nullptr ? search_
.refIndices_
[j
] : j
);
1036 const int refId
= search_
.refExclusionIds_
[index
];
1037 while (exclind_
< nexcl_
&& excl_
[exclind_
] < refId
)
1041 if (exclind_
< nexcl_
&& refId
== excl_
[exclind_
])
1050 void AnalysisNeighborhoodPairSearchImpl::startSearch(
1051 const AnalysisNeighborhoodPositions
&positions
)
1053 selfSearchMode_
= false;
1054 testPosCount_
= positions
.count_
;
1055 testPositions_
= positions
.x_
;
1056 testExclusionIds_
= positions
.exclusionIds_
;
1057 testIndices_
= positions
.indices_
;
1058 GMX_RELEASE_ASSERT(search_
.excls_
== nullptr || testExclusionIds_
!= nullptr,
1059 "Exclusion IDs must be set when exclusions are enabled");
1060 if (positions
.index_
< 0)
1066 // Somewhat of a hack: setup the array such that only the last position
1068 testPosCount_
= positions
.index_
+ 1;
1069 reset(positions
.index_
);
1073 void AnalysisNeighborhoodPairSearchImpl::startSelfSearch()
1075 selfSearchMode_
= true;
1076 testPosCount_
= search_
.nref_
;
1077 testPositions_
= search_
.xref_
;
1078 testExclusionIds_
= search_
.refExclusionIds_
;
1079 testIndices_
= search_
.refIndices_
;
1080 GMX_RELEASE_ASSERT(search_
.excls_
== nullptr || testIndices_
== nullptr,
1081 "Exclusion IDs not implemented with indexed ref positions");
1085 template <class Action
>
1086 bool AnalysisNeighborhoodPairSearchImpl::searchNext(Action action
)
1088 while (testIndex_
< testPosCount_
)
1092 int cai
= prevcai_
+ 1;
1097 const int ci
= search_
.shiftCell(currCell_
, shift
);
1098 if (selfSearchMode_
&& ci
> testCellIndex_
)
1102 const int cellSize
= static_cast<int>(search_
.cells_
[ci
].size());
1103 for (; cai
< cellSize
; ++cai
)
1105 const int i
= search_
.cells_
[ci
][cai
];
1106 if (selfSearchMode_
&& ci
== testCellIndex_
&& i
>= testIndex_
)
1115 rvec_sub(search_
.xref_
[i
], xtest_
, dx
);
1116 rvec_sub(dx
, shift
, dx
);
1119 ? dx
[XX
]*dx
[XX
] + dx
[YY
]*dx
[YY
]
1121 if (r2
<= search_
.cutoff2_
)
1123 if (action(i
, r2
, dx
))
1128 copy_rvec(dx
, prevdx_
);
1136 while (search_
.nextCell(testcell_
, currCell_
, cellBound_
));
1140 for (int i
= previ_
+ 1; i
< search_
.nref_
; ++i
)
1147 if (search_
.pbc_
.ePBC
!= epbcNONE
)
1149 pbc_dx(&search_
.pbc_
, search_
.xref_
[i
], xtest_
, dx
);
1153 rvec_sub(search_
.xref_
[i
], xtest_
, dx
);
1157 ? dx
[XX
]*dx
[XX
] + dx
[YY
]*dx
[YY
]
1159 if (r2
<= search_
.cutoff2_
)
1161 if (action(i
, r2
, dx
))
1165 copy_rvec(dx
, prevdx_
);
1176 void AnalysisNeighborhoodPairSearchImpl::initFoundPair(
1177 AnalysisNeighborhoodPair
*pair
) const
1181 *pair
= AnalysisNeighborhoodPair();
1185 *pair
= AnalysisNeighborhoodPair(previ_
, testIndex_
, prevr2_
, prevdx_
);
1189 } // namespace internal
1195 * Search action to find the next neighbor.
1197 * Used as the action for AnalysisNeighborhoodPairSearchImpl::searchNext() to
1198 * find the next neighbor.
1200 * Simply breaks the loop on the first found neighbor.
1202 bool withinAction(int /*i*/, real
/*r2*/, const rvec
/* dx */)
1208 * Search action find the minimum distance.
1210 * Used as the action for AnalysisNeighborhoodPairSearchImpl::searchNext() to
1211 * find the nearest neighbor.
1213 * With this action, AnalysisNeighborhoodPairSearchImpl::searchNext() always
1214 * returns false, and the output is put into the variables passed by pointer
1215 * into the constructor. If no neighbors are found, the output variables are
1216 * not modified, i.e., the caller must initialize them.
1222 * Initializes the action with given output locations.
1224 * \param[out] closestPoint Index of the closest reference location.
1225 * \param[out] minDist2 Minimum distance squared.
1226 * \param[out] dx Shortest distance vector.
1228 * The constructor call does not modify the pointed values, but only
1229 * stores the pointers for later use.
1230 * See the class description for additional semantics.
1232 MindistAction(int *closestPoint
, real
*minDist2
, rvec
*dx
)
1233 : closestPoint_(*closestPoint
), minDist2_(*minDist2
), dx_(*dx
)
1236 //! Copies the action.
1237 MindistAction(const MindistAction
&) = default;
1239 //! Processes a neighbor to find the nearest point.
1240 bool operator()(int i
, real r2
, const rvec dx
)
1256 GMX_DISALLOW_ASSIGN(MindistAction
);
1261 /********************************************************************
1262 * AnalysisNeighborhood::Impl
1265 class AnalysisNeighborhood::Impl
1268 typedef AnalysisNeighborhoodSearch::ImplPointer SearchImplPointer
;
1269 typedef std::vector
<SearchImplPointer
> SearchList
;
1272 : cutoff_(0), excls_(nullptr), mode_(eSearchMode_Automatic
), bXY_(false)
1277 SearchList::const_iterator i
;
1278 for (i
= searchList_
.begin(); i
!= searchList_
.end(); ++i
)
1280 GMX_RELEASE_ASSERT(i
->unique(),
1281 "Dangling AnalysisNeighborhoodSearch reference");
1285 SearchImplPointer
getSearch();
1287 Mutex createSearchMutex_
;
1288 SearchList searchList_
;
1290 const t_blocka
*excls_
;
1295 AnalysisNeighborhood::Impl::SearchImplPointer
1296 AnalysisNeighborhood::Impl::getSearch()
1298 lock_guard
<Mutex
> lock(createSearchMutex_
);
1299 // TODO: Consider whether this needs to/can be faster, e.g., by keeping a
1300 // separate pool of unused search objects.
1301 SearchList::const_iterator i
;
1302 for (i
= searchList_
.begin(); i
!= searchList_
.end(); ++i
)
1309 SearchImplPointer
search(new internal::AnalysisNeighborhoodSearchImpl(cutoff_
));
1310 searchList_
.push_back(search
);
1314 /********************************************************************
1315 * AnalysisNeighborhood
1318 AnalysisNeighborhood::AnalysisNeighborhood()
1323 AnalysisNeighborhood::~AnalysisNeighborhood()
1327 void AnalysisNeighborhood::setCutoff(real cutoff
)
1329 GMX_RELEASE_ASSERT(impl_
->searchList_
.empty(),
1330 "Changing the cutoff after initSearch() not currently supported");
1331 impl_
->cutoff_
= cutoff
;
1334 void AnalysisNeighborhood::setXYMode(bool bXY
)
1339 void AnalysisNeighborhood::setTopologyExclusions(const t_blocka
*excls
)
1341 GMX_RELEASE_ASSERT(impl_
->searchList_
.empty(),
1342 "Changing the exclusions after initSearch() not currently supported");
1343 impl_
->excls_
= excls
;
1346 void AnalysisNeighborhood::setMode(SearchMode mode
)
1348 impl_
->mode_
= mode
;
1351 AnalysisNeighborhood::SearchMode
AnalysisNeighborhood::mode() const
1353 return impl_
->mode_
;
1356 AnalysisNeighborhoodSearch
1357 AnalysisNeighborhood::initSearch(const t_pbc
*pbc
,
1358 const AnalysisNeighborhoodPositions
&positions
)
1360 Impl::SearchImplPointer
search(impl_
->getSearch());
1361 search
->init(mode(), impl_
->bXY_
, impl_
->excls_
,
1363 return AnalysisNeighborhoodSearch(search
);
1366 /********************************************************************
1367 * AnalysisNeighborhoodSearch
1370 AnalysisNeighborhoodSearch::AnalysisNeighborhoodSearch()
1374 AnalysisNeighborhoodSearch::AnalysisNeighborhoodSearch(const ImplPointer
&impl
)
1379 void AnalysisNeighborhoodSearch::reset()
1384 AnalysisNeighborhood::SearchMode
AnalysisNeighborhoodSearch::mode() const
1386 GMX_RELEASE_ASSERT(impl_
, "Accessing an invalid search object");
1387 return (impl_
->usesGridSearch()
1388 ? AnalysisNeighborhood::eSearchMode_Grid
1389 : AnalysisNeighborhood::eSearchMode_Simple
);
1392 bool AnalysisNeighborhoodSearch::isWithin(
1393 const AnalysisNeighborhoodPositions
&positions
) const
1395 GMX_RELEASE_ASSERT(impl_
, "Accessing an invalid search object");
1396 internal::AnalysisNeighborhoodPairSearchImpl
pairSearch(*impl_
);
1397 pairSearch
.startSearch(positions
);
1398 return pairSearch
.searchNext(&withinAction
);
1401 real
AnalysisNeighborhoodSearch::minimumDistance(
1402 const AnalysisNeighborhoodPositions
&positions
) const
1404 GMX_RELEASE_ASSERT(impl_
, "Accessing an invalid search object");
1405 internal::AnalysisNeighborhoodPairSearchImpl
pairSearch(*impl_
);
1406 pairSearch
.startSearch(positions
);
1407 real minDist2
= impl_
->cutoffSquared();
1408 int closestPoint
= -1;
1409 rvec dx
= {0.0, 0.0, 0.0};
1410 MindistAction
action(&closestPoint
, &minDist2
, &dx
);
1411 (void)pairSearch
.searchNext(action
);
1412 return std::sqrt(minDist2
);
1415 AnalysisNeighborhoodPair
1416 AnalysisNeighborhoodSearch::nearestPoint(
1417 const AnalysisNeighborhoodPositions
&positions
) const
1419 GMX_RELEASE_ASSERT(impl_
, "Accessing an invalid search object");
1420 internal::AnalysisNeighborhoodPairSearchImpl
pairSearch(*impl_
);
1421 pairSearch
.startSearch(positions
);
1422 real minDist2
= impl_
->cutoffSquared();
1423 int closestPoint
= -1;
1424 rvec dx
= {0.0, 0.0, 0.0};
1425 MindistAction
action(&closestPoint
, &minDist2
, &dx
);
1426 (void)pairSearch
.searchNext(action
);
1427 return AnalysisNeighborhoodPair(closestPoint
, 0, minDist2
, dx
);
1430 AnalysisNeighborhoodPairSearch
1431 AnalysisNeighborhoodSearch::startSelfPairSearch() const
1433 GMX_RELEASE_ASSERT(impl_
, "Accessing an invalid search object");
1434 Impl::PairSearchImplPointer
pairSearch(impl_
->getPairSearch());
1435 pairSearch
->startSelfSearch();
1436 return AnalysisNeighborhoodPairSearch(pairSearch
);
1439 AnalysisNeighborhoodPairSearch
1440 AnalysisNeighborhoodSearch::startPairSearch(
1441 const AnalysisNeighborhoodPositions
&positions
) const
1443 GMX_RELEASE_ASSERT(impl_
, "Accessing an invalid search object");
1444 Impl::PairSearchImplPointer
pairSearch(impl_
->getPairSearch());
1445 pairSearch
->startSearch(positions
);
1446 return AnalysisNeighborhoodPairSearch(pairSearch
);
1449 /********************************************************************
1450 * AnalysisNeighborhoodPairSearch
1453 AnalysisNeighborhoodPairSearch::AnalysisNeighborhoodPairSearch(
1454 const ImplPointer
&impl
)
1459 bool AnalysisNeighborhoodPairSearch::findNextPair(AnalysisNeighborhoodPair
*pair
)
1461 bool bFound
= impl_
->searchNext(&withinAction
);
1462 impl_
->initFoundPair(pair
);
1466 void AnalysisNeighborhoodPairSearch::skipRemainingPairsForTestPosition()
1468 impl_
->nextTestPosition();