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 * - 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
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"
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
)
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
])
103 if (maxBound
[d
] < x
[i
][d
])
105 maxBound
[d
] = x
[i
][d
];
109 rvec_sub(maxBound
, origin
, size
);
117 /********************************************************************
118 * Implementation class declarations
121 class AnalysisNeighborhoodSearchImpl
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();
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 t_blocka
*excls
,
145 const AnalysisNeighborhoodPositions
&positions
);
146 PairSearchImplPointer
getPairSearch();
148 real
cutoffSquared() const { return cutoff2_
; }
149 bool usesGridSearch() const { return bGrid_
; }
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
);
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
167 * \returns `false` if grid search is not suitable.
169 bool initGridCells(const matrix box
, bool bSingleCell
[DIM
],
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
);
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;
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;
202 * Adds an index into a grid cell.
204 * \param[in] cell Fractional cell coordinates into which \p i should
206 * \param[in] i Index to add.
208 * \p cell should satisfy the conditions that \p mapPointToGridCell()
211 void addToGridCell(const rvec cell
, int i
);
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
234 void initCellRange(const rvec centerCell
, ivec cell
,
235 ivec upperBound
, int dim
) const;
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;
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.
259 //! The cutoff squared.
261 //! Whether to do searching in XY plane only.
264 //! Number of reference points for the current frame.
266 //! Reference point positions.
268 //! Reference position exclusion IDs.
269 const int *refExclusionIds_
;
270 //! Reference position indices (NULL if no indices).
271 const int *refIndices_
;
273 const t_blocka
*excls_
;
277 //! Whether grid searching is actually used for the current positions.
279 //! false if the box is rectangular.
281 //! Whether the grid is periodic in a dimension.
283 //! Array for storing in-unit-cell reference positions.
284 std::vector
<RVec
> xrefAlloc_
;
285 //! Origin of the grid (zero for periodic dimensions).
287 //! Size of a single grid cell.
289 //! Inverse of \p cellSize_. Zero for dimensions where grid is not used.
292 * Shift in cell coordinates (for triclinic boxes) in X when crossing
293 * the Z periodic boundary.
297 * Shift in cell coordinates (for triclinic boxes) in Y when crossing
298 * the Z periodic boundary.
302 * Shift in cell coordinates (for triclinic boxes) in X when crossing
303 * the Y periodic boundary.
306 //! Number of cells along each dimension.
308 //! Data structure to hold the grid cell contents.
311 Mutex createPairSearchMutex_
;
312 PairSearchList pairSearchList_
;
314 friend class AnalysisNeighborhoodPairSearchImpl
;
316 GMX_DISALLOW_COPY_AND_ASSIGN(AnalysisNeighborhoodSearchImpl
);
319 class AnalysisNeighborhoodPairSearchImpl
322 explicit AnalysisNeighborhoodPairSearchImpl(const AnalysisNeighborhoodSearchImpl
&search
)
326 testPositions_
= nullptr;
327 testExclusionIds_
= nullptr;
328 testIndices_
= nullptr;
332 clear_rvec(testcell_
);
333 clear_ivec(currCell_
);
334 clear_ivec(cellBound_
);
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();
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.
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.
366 //! Exclusions for current test particle.
368 //! Index of the currently active test position in \p testPositions_.
370 //! Stores test position during a pair loop.
372 //! Stores the previous returned position during a pair loop.
374 //! Stores the pair distance corresponding to previ_;
376 //! Stores the shortest distance vector corresponding to previ_;
378 //! Stores the current exclusion index during loops.
380 //! Stores the fractional test particle cell location during loops.
382 //! Stores the current cell during pair loops.
384 //! Stores the current loop upper bounds for each dimension during pair loops.
386 //! Stores the index within the current cell during pair loops.
389 GMX_DISALLOW_COPY_AND_ASSIGN(AnalysisNeighborhoodPairSearchImpl
);
392 /********************************************************************
393 * AnalysisNeighborhoodSearchImpl
396 AnalysisNeighborhoodSearchImpl::AnalysisNeighborhoodSearchImpl(real cutoff
)
402 cutoff_
= cutoff2_
= GMX_REAL_MAX
;
407 cutoff2_
= gmx::square(cutoff_
);
412 refExclusionIds_
= nullptr;
413 refIndices_
= nullptr;
414 std::memset(&pbc_
, 0, sizeof(pbc_
));
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
)
452 PairSearchImplPointer
pairSearch(new AnalysisNeighborhoodPairSearchImpl(*this));
453 pairSearchList_
.push_back(pairSearch
);
457 bool AnalysisNeighborhoodSearchImpl::checkGridSearchEfficiency(bool bForce
)
459 // Find the extent of the sphere in cells.
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;
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
482 coveredCells
*= coveredCount
;
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
;
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
)
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;
532 for (int dd
= 0; dd
< DIM
; ++dd
)
534 const real boxSize
= box
[dd
][dd
];
535 if (boxSize
< targetsize
)
537 bSingleCell
[dd
] = true;
552 if (dimCount
== 0 || dimCount
== prevDimCount
)
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
)
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)
579 totalCellCount
*= cellCount
;
580 ncelldim_
[dd
] = cellCount
;
582 if (totalCellCount
<= 3)
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
)
600 bool AnalysisNeighborhoodSearchImpl::initGrid(
601 const t_pbc
&pbc
, int posCount
, const rvec x
[], bool bForce
)
611 bGridPBC_
[XX
] = false;
612 bGridPBC_
[YY
] = false;
613 bGridPBC_
[ZZ
] = false;
616 bGridPBC_
[XX
] = true;
617 bGridPBC_
[YY
] = true;
618 bGridPBC_
[ZZ
] = false;
621 bGridPBC_
[XX
] = true;
622 bGridPBC_
[YY
] = true;
623 bGridPBC_
[ZZ
] = true;
626 // Grid searching not supported for now with screw.
630 bool bSingleCell
[DIM
] = {false, false, bXY_
};
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
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
];
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;
658 if (!initGridCells(box
, bSingleCell
, posCount
))
663 bTric_
= TRICLINIC(pbc
.box
);
664 for (int dd
= 0; dd
< DIM
; ++dd
)
666 cellSize_
[dd
] = box
[dd
][dd
] / ncelldim_
[dd
];
669 invCellSize_
[dd
] = 0.0;
673 invCellSize_
[dd
] = 1.0 / cellSize_
[dd
];
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
,
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
];
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
)
731 for (int dd
= 0; dd
< DIM
; ++dd
)
733 int cellIndex
= static_cast<int>(floor(cell
[dd
]));
736 const int cellCount
= ncelldim_
[dd
];
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
;
766 if (currCell
[ZZ
] < 0)
768 startOffset
+= cellShiftZY_
;
769 endOffset
+= cellShiftZY_
;
771 else if (currCell
[ZZ
] >= ncelldim_
[ZZ
])
773 startOffset
-= cellShiftZY_
;
774 endOffset
-= cellShiftZY_
;
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_
;
801 // For non-periodic dimensions, clamp to the actual grid edges.
804 // If endOffset < 0 or startOffset > N, these may cause the whole
805 // test position/grid plane/grid row to be skipped.
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
828 if (cell
[dim
] > upperBound
[dim
])
833 for (int d
= dim
- 1; d
>= 0; --d
)
835 initCellRange(centerCell
, cell
, upperBound
, d
);
836 if (cell
[d
] > upperBound
[d
])
847 int AnalysisNeighborhoodSearchImpl::shiftCell(const ivec cell
, rvec shift
) const
850 copy_ivec(cell
, shiftedCell
);
853 for (int d
= 0; d
< DIM
; ++d
)
855 const int cellCount
= ncelldim_
[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
,
880 const t_blocka
*excls
,
882 const AnalysisNeighborhoodPositions
&positions
)
884 GMX_RELEASE_ASSERT(positions
.index_
== -1,
885 "Individual indexed positions not supported as reference");
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.
904 copy_mat(pbc
->box
, box
);
906 set_pbc(&pbc_
, epbcXY
, box
);
908 else if (pbc
!= nullptr)
914 pbc_
.ePBC
= epbcNONE
;
917 nref_
= positions
.count_
;
918 if (mode
== AnalysisNeighborhood::eSearchMode_Simple
)
924 bGrid_
= initGrid(pbc_
, positions
.count_
, positions
.x_
,
925 mode
== AnalysisNeighborhood::eSearchMode_Grid
);
927 refIndices_
= positions
.indices_
;
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
;
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
]);
952 xref_
= positions
.x_
;
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_
)
976 (testIndices_
!= nullptr ? testIndices_
[testIndex
] : testIndex
);
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
);
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
];
1006 clear_rvec(prevdx_
);
1011 void AnalysisNeighborhoodPairSearchImpl::nextTestPosition()
1013 if (testIndex_
< testPosCount_
)
1020 bool AnalysisNeighborhoodPairSearchImpl::isExcluded(int j
)
1022 if (exclind_
< nexcl_
)
1025 (search_
.refIndices_
!= nullptr ? search_
.refIndices_
[j
] : j
);
1026 const int refId
= search_
.refExclusionIds_
[index
];
1027 while (exclind_
< nexcl_
&& excl_
[exclind_
] < refId
)
1031 if (exclind_
< nexcl_
&& refId
== excl_
[exclind_
])
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)
1055 // Somewhat of a hack: setup the array such that only the last position
1057 testPosCount_
= positions
.index_
+ 1;
1058 reset(positions
.index_
);
1062 template <class Action
>
1063 bool AnalysisNeighborhoodPairSearchImpl::searchNext(Action action
)
1065 while (testIndex_
< testPosCount_
)
1069 int cai
= prevcai_
+ 1;
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
];
1084 rvec_sub(search_
.xref_
[i
], xtest_
, dx
);
1085 rvec_sub(dx
, shift
, dx
);
1088 ? dx
[XX
]*dx
[XX
] + dx
[YY
]*dx
[YY
]
1090 if (r2
<= search_
.cutoff2_
)
1092 if (action(i
, r2
, dx
))
1097 copy_rvec(dx
, prevdx_
);
1105 while (search_
.nextCell(testcell_
, currCell_
, cellBound_
));
1109 for (int i
= previ_
+ 1; i
< search_
.nref_
; ++i
)
1116 if (search_
.pbc_
.ePBC
!= epbcNONE
)
1118 pbc_dx(&search_
.pbc_
, search_
.xref_
[i
], xtest_
, dx
);
1122 rvec_sub(search_
.xref_
[i
], xtest_
, dx
);
1126 ? dx
[XX
]*dx
[XX
] + dx
[YY
]*dx
[YY
]
1128 if (r2
<= search_
.cutoff2_
)
1130 if (action(i
, r2
, dx
))
1134 copy_rvec(dx
, prevdx_
);
1145 void AnalysisNeighborhoodPairSearchImpl::initFoundPair(
1146 AnalysisNeighborhoodPair
*pair
) const
1150 *pair
= AnalysisNeighborhoodPair();
1154 *pair
= AnalysisNeighborhoodPair(previ_
, testIndex_
, prevr2_
, prevdx_
);
1158 } // namespace internal
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 */)
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.
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
)
1225 GMX_DISALLOW_ASSIGN(MindistAction
);
1230 /********************************************************************
1231 * AnalysisNeighborhood::Impl
1234 class AnalysisNeighborhood::Impl
1237 typedef AnalysisNeighborhoodSearch::ImplPointer SearchImplPointer
;
1238 typedef std::vector
<SearchImplPointer
> SearchList
;
1241 : cutoff_(0), excls_(nullptr), mode_(eSearchMode_Automatic
), bXY_(false)
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_
;
1259 const t_blocka
*excls_
;
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
)
1278 SearchImplPointer
search(new internal::AnalysisNeighborhoodSearchImpl(cutoff_
));
1279 searchList_
.push_back(search
);
1283 /********************************************************************
1284 * AnalysisNeighborhood
1287 AnalysisNeighborhood::AnalysisNeighborhood()
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
)
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_
,
1332 return AnalysisNeighborhoodSearch(search
);
1335 /********************************************************************
1336 * AnalysisNeighborhoodSearch
1339 AnalysisNeighborhoodSearch::AnalysisNeighborhoodSearch()
1343 AnalysisNeighborhoodSearch::AnalysisNeighborhoodSearch(const ImplPointer
&impl
)
1348 void AnalysisNeighborhoodSearch::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
)
1419 bool AnalysisNeighborhoodPairSearch::findNextPair(AnalysisNeighborhoodPair
*pair
)
1421 bool bFound
= impl_
->searchNext(&withinAction
);
1422 impl_
->initFoundPair(pair
);
1426 void AnalysisNeighborhoodPairSearch::skipRemainingPairsForTestPosition()
1428 impl_
->nextTestPosition();