2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2009,2010,2011,2012,2013,2014, 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).
40 * The grid implementation could still be optimized in several different ways:
41 * - Triclinic grid cells are not the most efficient shape, but make PBC
43 * - Precalculating the required PBC shift for a pair of cells outside the
44 * inner loop. After this is done, it should be quite straightforward to
45 * move to rectangular cells.
46 * - Pruning grid cells from the search list if they are completely outside
47 * the sphere that is being considered.
48 * - A better heuristic could be added for falling back to simple loops for a
49 * small number of reference particles.
50 * - A better heuristic for selecting the grid size.
51 * - A multi-level grid implementation could be used to be able to use small
52 * grids for short cutoffs with very inhomogeneous particle distributions
53 * without a memory cost.
55 * \author Teemu Murtola <teemu.murtola@gmail.com>
56 * \ingroup module_selection
58 #include "gromacs/selection/nbsearch.h"
66 #include "thread_mpi/mutex.h"
68 #include "gromacs/legacyheaders/names.h"
70 #include "gromacs/math/vec.h"
71 #include "gromacs/pbcutil/pbc.h"
72 #include "gromacs/selection/position.h"
73 #include "gromacs/topology/block.h"
74 #include "gromacs/utility/arrayref.h"
75 #include "gromacs/utility/exceptions.h"
76 #include "gromacs/utility/gmxassert.h"
77 #include "gromacs/utility/smalloc.h"
78 #include "gromacs/utility/stringutil.h"
86 /********************************************************************
87 * Implementation class declarations
90 class AnalysisNeighborhoodSearchImpl
93 typedef AnalysisNeighborhoodPairSearch::ImplPointer
94 PairSearchImplPointer;
95 typedef std::vector<PairSearchImplPointer> PairSearchList;
96 typedef std::vector<std::vector<int> > CellList;
98 explicit AnalysisNeighborhoodSearchImpl(real cutoff);
99 ~AnalysisNeighborhoodSearchImpl();
102 * Initializes the search with a given box and reference positions.
104 * \param[in] mode Search mode to use.
105 * \param[in] bXY Whether to use 2D searching.
106 * \param[in] excls Exclusions.
107 * \param[in] pbc PBC information.
108 * \param[in] positions Set of reference positions.
110 void init(AnalysisNeighborhood::SearchMode mode,
112 const t_blocka *excls,
114 const AnalysisNeighborhoodPositions &positions);
115 PairSearchImplPointer getPairSearch();
117 real cutoffSquared() const { return cutoff2_; }
118 bool usesGridSearch() const { return bGrid_; }
121 //! Calculates offsets to neighboring grid cells that should be considered.
122 void initGridCellNeighborList();
124 * Determines a suitable grid size and sets up the cells.
126 * \param[in] pbc Information about the box.
127 * \returns false if grid search is not suitable.
129 bool initGridCells(const t_pbc &pbc);
131 * Sets ua a search grid for a given box.
133 * \param[in] pbc Information about the box.
134 * \returns false if grid search is not suitable.
136 bool initGrid(const t_pbc &pbc);
138 * Maps a point into a grid cell.
140 * \param[in] x Point to map.
141 * \param[out] cell Indices of the grid cell in which \p x lies.
143 * \p x should be within the triclinic unit cell.
145 void mapPointToGridCell(const rvec x, ivec cell) const;
147 * Calculates linear index of a grid cell.
149 * \param[in] cell Cell indices.
150 * \returns Linear index of \p cell.
152 int getGridCellIndex(const ivec cell) const;
154 * Adds an index into a grid cell.
156 * \param[in] cell Cell into which \p i should be added.
157 * \param[in] i Index to add.
159 void addToGridCell(const ivec cell, int i);
161 //! Whether to try grid searching.
165 //! The cutoff squared.
167 //! Whether to do searching in XY plane only.
170 //! Number of reference points for the current frame.
172 //! Reference point positions.
174 //! Reference position exclusion IDs.
175 const int *refExclusionIds_;
177 const t_blocka *excls_;
181 //! Whether grid searching is actually used for the current positions.
183 //! Array allocated for storing in-unit-cell reference positions.
185 //! Allocation count for xref_alloc.
187 //! false if the box is rectangular.
189 //! Box vectors of a single grid cell.
191 //! The reciprocal cell vectors as columns; the inverse of \p cellbox.
193 //! Number of cells along each dimension.
195 //! Data structure to hold the grid cell contents.
197 //! Number of neighboring cells to consider.
199 //! Offsets of the neighboring cells to consider.
201 //! Allocation count for \p gnboffs.
204 tMPI::mutex createPairSearchMutex_;
205 PairSearchList pairSearchList_;
207 friend class AnalysisNeighborhoodPairSearchImpl;
209 GMX_DISALLOW_COPY_AND_ASSIGN(AnalysisNeighborhoodSearchImpl);
212 class AnalysisNeighborhoodPairSearchImpl
215 explicit AnalysisNeighborhoodPairSearchImpl(const AnalysisNeighborhoodSearchImpl &search)
218 testExclusionIds_ = NULL;
222 clear_ivec(testcell_);
226 //! Initializes a search to find reference positions neighboring \p x.
227 void startSearch(const AnalysisNeighborhoodPositions &positions);
228 //! Searches for the next neighbor.
229 template <class Action>
230 bool searchNext(Action action);
231 //! Initializes a pair representing the pair found by searchNext().
232 void initFoundPair(AnalysisNeighborhoodPair *pair) const;
233 //! Advances to the next test position, skipping any remaining pairs.
234 void nextTestPosition();
237 //! Clears the loop indices.
238 void reset(int testIndex);
239 //! Checks whether a reference positiong should be excluded.
240 bool isExcluded(int j);
242 //! Parent search object.
243 const AnalysisNeighborhoodSearchImpl &search_;
244 //! Reference to the test positions.
245 ConstArrayRef<rvec> testPositions_;
246 //! Reference to the test exclusion indices.
247 const int *testExclusionIds_;
248 //! Number of excluded reference positions for current test particle.
250 //! Exclusions for current test particle.
252 //! Index of the currently active test position in \p testPositions_.
254 //! Stores test position during a pair loop.
256 //! Stores the previous returned position during a pair loop.
258 //! Stores the pair distance corresponding to previ_;
260 //! Stores the current exclusion index during loops.
262 //! Stores the test particle cell index during loops.
264 //! Stores the current cell neighbor index during pair loops.
266 //! Stores the index within the current cell during pair loops.
269 GMX_DISALLOW_COPY_AND_ASSIGN(AnalysisNeighborhoodPairSearchImpl);
272 /********************************************************************
273 * AnalysisNeighborhoodSearchImpl
276 AnalysisNeighborhoodSearchImpl::AnalysisNeighborhoodSearchImpl(real cutoff)
282 cutoff_ = GMX_REAL_MAX;
285 cutoff2_ = sqr(cutoff_);
290 refExclusionIds_ = NULL;
291 std::memset(&pbc_, 0, sizeof(pbc_));
299 clear_mat(recipcell_);
300 clear_ivec(ncelldim_);
307 AnalysisNeighborhoodSearchImpl::~AnalysisNeighborhoodSearchImpl()
309 PairSearchList::const_iterator i;
310 for (i = pairSearchList_.begin(); i != pairSearchList_.end(); ++i)
312 GMX_RELEASE_ASSERT(i->unique(),
313 "Dangling AnalysisNeighborhoodPairSearch reference");
319 AnalysisNeighborhoodSearchImpl::PairSearchImplPointer
320 AnalysisNeighborhoodSearchImpl::getPairSearch()
322 tMPI::lock_guard<tMPI::mutex> lock(createPairSearchMutex_);
323 // TODO: Consider whether this needs to/can be faster, e.g., by keeping a
324 // separate pool of unused search objects.
325 PairSearchList::const_iterator i;
326 for (i = pairSearchList_.begin(); i != pairSearchList_.end(); ++i)
333 PairSearchImplPointer pairSearch(new AnalysisNeighborhoodPairSearchImpl(*this));
334 pairSearchList_.push_back(pairSearch);
338 void AnalysisNeighborhoodSearchImpl::initGridCellNeighborList()
340 int maxx, maxy, maxz;
343 /* Find the extent of the sphere in triclinic coordinates */
344 maxz = static_cast<int>(cutoff_ * recipcell_[ZZ][ZZ]) + 1;
345 rvnorm = sqrt(sqr(recipcell_[YY][YY]) + sqr(recipcell_[ZZ][YY]));
346 maxy = static_cast<int>(cutoff_ * rvnorm) + 1;
347 rvnorm = sqrt(sqr(recipcell_[XX][XX]) + sqr(recipcell_[YY][XX])
348 + sqr(recipcell_[ZZ][XX]));
349 maxx = static_cast<int>(cutoff_ * rvnorm) + 1;
351 /* Calculate the number of cells and reallocate if necessary */
352 ngridnb_ = (2 * maxx + 1) * (2 * maxy + 1) * (2 * maxz + 1);
353 if (gnboffs_nalloc_ < ngridnb_)
355 gnboffs_nalloc_ = ngridnb_;
356 srenew(gnboffs_, gnboffs_nalloc_);
359 /* Store the whole cube */
360 /* TODO: Prune off corners that are not needed */
362 for (int x = -maxx; x <= maxx; ++x)
364 for (int y = -maxy; y <= maxy; ++y)
366 for (int z = -maxz; z <= maxz; ++z)
377 bool AnalysisNeighborhoodSearchImpl::initGridCells(const t_pbc &pbc)
379 const real targetsize =
380 pow(pbc.box[XX][XX] * pbc.box[YY][YY] * pbc.box[ZZ][ZZ]
381 * 10 / nref_, static_cast<real>(1./3.));
384 for (int dd = 0; dd < DIM; ++dd)
386 ncelldim_[dd] = static_cast<int>(pbc.box[dd][dd] / targetsize);
387 cellCount *= ncelldim_[dd];
388 if (ncelldim_[dd] < 3)
393 // Never decrease the size of the cell vector to avoid reallocating
394 // memory for the nested vectors. The actual size of the vector is not
395 // used outside this function.
396 if (cells_.size() < static_cast<size_t>(cellCount))
398 cells_.resize(cellCount);
400 for (int ci = 0; ci < cellCount; ++ci)
407 bool AnalysisNeighborhoodSearchImpl::initGrid(const t_pbc &pbc)
409 /* TODO: This check could be improved. */
410 if (0.5*pbc.max_cutoff2 < cutoff2_)
415 if (!initGridCells(pbc))
420 bTric_ = TRICLINIC(pbc.box);
423 for (int dd = 0; dd < DIM; ++dd)
425 svmul(1.0 / ncelldim_[dd], pbc.box[dd], cellbox_[dd]);
427 m_inv_ur0(cellbox_, recipcell_);
431 for (int dd = 0; dd < DIM; ++dd)
433 cellbox_[dd][dd] = pbc.box[dd][dd] / ncelldim_[dd];
434 recipcell_[dd][dd] = 1.0 / cellbox_[dd][dd];
437 initGridCellNeighborList();
441 void AnalysisNeighborhoodSearchImpl::mapPointToGridCell(const rvec x,
448 tmvmul_ur0(recipcell_, x, tx);
449 for (int dd = 0; dd < DIM; ++dd)
451 cell[dd] = static_cast<int>(tx[dd]);
456 for (int dd = 0; dd < DIM; ++dd)
458 cell[dd] = static_cast<int>(x[dd] * recipcell_[dd][dd]);
463 int AnalysisNeighborhoodSearchImpl::getGridCellIndex(const ivec cell) const
465 GMX_ASSERT(cell[XX] >= 0 && cell[XX] < ncelldim_[XX],
466 "Grid cell X index out of range");
467 GMX_ASSERT(cell[YY] >= 0 && cell[YY] < ncelldim_[YY],
468 "Grid cell Y index out of range");
469 GMX_ASSERT(cell[ZZ] >= 0 && cell[ZZ] < ncelldim_[ZZ],
470 "Grid cell Z index out of range");
471 return cell[XX] + cell[YY] * ncelldim_[XX]
472 + cell[ZZ] * ncelldim_[XX] * ncelldim_[YY];
475 void AnalysisNeighborhoodSearchImpl::addToGridCell(const ivec cell, int i)
477 const int ci = getGridCellIndex(cell);
478 cells_[ci].push_back(i);
481 void AnalysisNeighborhoodSearchImpl::init(
482 AnalysisNeighborhood::SearchMode mode,
484 const t_blocka *excls,
486 const AnalysisNeighborhoodPositions &positions)
488 GMX_RELEASE_ASSERT(positions.index_ == -1,
489 "Individual indexed positions not supported as reference");
491 if (bXY_ && pbc->ePBC != epbcNONE)
493 if (pbc->ePBC != epbcXY && pbc->ePBC != epbcXYZ)
495 std::string message =
496 formatString("Computations in the XY plane are not supported with PBC type '%s'",
498 GMX_THROW(NotImplementedError(message));
500 if (std::fabs(pbc->box[ZZ][XX]) > GMX_REAL_EPS*pbc->box[ZZ][ZZ] ||
501 std::fabs(pbc->box[ZZ][YY]) > GMX_REAL_EPS*pbc->box[ZZ][ZZ])
503 GMX_THROW(NotImplementedError("Computations in the XY plane are not supported when the last box vector is not parallel to the Z axis"));
505 set_pbc(&pbc_, epbcXY, const_cast<rvec *>(pbc->box));
507 else if (pbc != NULL)
513 pbc_.ePBC = epbcNONE;
515 nref_ = positions.count_;
516 // TODO: Consider whether it would be possible to support grid searching in
518 if (mode == AnalysisNeighborhood::eSearchMode_Simple
519 || pbc_.ePBC != epbcXYZ)
525 // TODO: Actually implement forcing eSearchMode_Grid
526 bGrid_ = initGrid(pbc_);
530 if (xref_nalloc_ < nref_)
532 srenew(xref_alloc_, nref_);
533 xref_nalloc_ = nref_;
537 for (int i = 0; i < nref_; ++i)
539 copy_rvec(positions.x_[i], xref_alloc_[i]);
541 put_atoms_in_triclinic_unitcell(ecenterTRIC, pbc_.box,
543 for (int i = 0; i < nref_; ++i)
547 mapPointToGridCell(xref_[i], refcell);
548 addToGridCell(refcell, i);
553 xref_ = positions.x_;
556 refExclusionIds_ = NULL;
559 // TODO: Check that the IDs are ascending, or remove the limitation.
560 refExclusionIds_ = positions.exclusionIds_;
561 GMX_RELEASE_ASSERT(refExclusionIds_ != NULL,
562 "Exclusion IDs must be set for reference positions "
563 "when exclusions are enabled");
567 /********************************************************************
568 * AnalysisNeighborhoodPairSearchImpl
571 void AnalysisNeighborhoodPairSearchImpl::reset(int testIndex)
573 testIndex_ = testIndex;
574 if (testIndex_ >= 0 && testIndex_ < static_cast<int>(testPositions_.size()))
576 copy_rvec(testPositions_[testIndex_], xtest_);
579 put_atoms_in_triclinic_unitcell(ecenterTRIC,
580 const_cast<rvec *>(search_.pbc_.box),
582 search_.mapPointToGridCell(xtest_, testcell_);
584 if (search_.excls_ != NULL)
586 const int exclIndex = testExclusionIds_[testIndex];
587 if (exclIndex < search_.excls_->nr)
589 const int startIndex = search_.excls_->index[exclIndex];
590 nexcl_ = search_.excls_->index[exclIndex + 1] - startIndex;
591 excl_ = &search_.excls_->a[startIndex];
607 void AnalysisNeighborhoodPairSearchImpl::nextTestPosition()
609 if (testIndex_ < static_cast<int>(testPositions_.size()))
616 bool AnalysisNeighborhoodPairSearchImpl::isExcluded(int j)
618 if (exclind_ < nexcl_)
620 const int refId = search_.refExclusionIds_[j];
621 while (exclind_ < nexcl_ && excl_[exclind_] < refId)
625 if (exclind_ < nexcl_ && refId == excl_[exclind_])
634 void AnalysisNeighborhoodPairSearchImpl::startSearch(
635 const AnalysisNeighborhoodPositions &positions)
637 testExclusionIds_ = positions.exclusionIds_;
638 GMX_RELEASE_ASSERT(search_.excls_ == NULL || testExclusionIds_ != NULL,
639 "Exclusion IDs must be set when exclusions are enabled");
640 if (positions.index_ < 0)
642 testPositions_ = ConstArrayRef<rvec>(positions.x_, positions.count_);
647 // Somewhat of a hack: setup the array such that only the last position
649 testPositions_ = ConstArrayRef<rvec>(positions.x_, positions.index_ + 1);
650 reset(positions.index_);
654 template <class Action>
655 bool AnalysisNeighborhoodPairSearchImpl::searchNext(Action action)
657 while (testIndex_ < static_cast<int>(testPositions_.size()))
661 GMX_RELEASE_ASSERT(!search_.bXY_, "Grid-based XY searches not implemented");
664 int cai = prevcai_ + 1;
666 for (; nbi < search_.ngridnb_; ++nbi)
670 ivec_add(testcell_, search_.gnboffs_[nbi], cell);
671 cell[XX] = (cell[XX] + search_.ncelldim_[XX]) % search_.ncelldim_[XX];
672 cell[YY] = (cell[YY] + search_.ncelldim_[YY]) % search_.ncelldim_[YY];
673 cell[ZZ] = (cell[ZZ] + search_.ncelldim_[ZZ]) % search_.ncelldim_[ZZ];
675 const int ci = search_.getGridCellIndex(cell);
676 const int cellSize = static_cast<int>(search_.cells_[ci].size());
677 /* TODO: Calculate the required PBC shift outside the inner loop */
678 for (; cai < cellSize; ++cai)
680 const int i = search_.cells_[ci][cai];
686 pbc_dx_aiuc(&search_.pbc_, xtest_, search_.xref_[i], dx);
687 const real r2 = norm2(dx);
688 if (r2 <= search_.cutoff2_)
706 for (int i = previ_ + 1; i < search_.nref_; ++i)
713 if (search_.pbc_.ePBC != epbcNONE)
715 pbc_dx(&search_.pbc_, xtest_, search_.xref_[i], dx);
719 rvec_sub(xtest_, search_.xref_[i], dx);
723 ? dx[XX]*dx[XX] + dx[YY]*dx[YY]
725 if (r2 <= search_.cutoff2_)
741 void AnalysisNeighborhoodPairSearchImpl::initFoundPair(
742 AnalysisNeighborhoodPair *pair) const
746 *pair = AnalysisNeighborhoodPair();
750 *pair = AnalysisNeighborhoodPair(previ_, testIndex_, prevr2_);
754 } // namespace internal
760 * Search action to find the next neighbor.
762 * Used as the action for AnalysisNeighborhoodPairSearchImpl::searchNext() to
763 * find the next neighbor.
765 * Simply breaks the loop on the first found neighbor.
767 bool withinAction(int /*i*/, real /*r2*/)
773 * Search action find the minimum distance.
775 * Used as the action for AnalysisNeighborhoodPairSearchImpl::searchNext() to
776 * find the nearest neighbor.
778 * With this action, AnalysisNeighborhoodPairSearchImpl::searchNext() always
779 * returns false, and the output is put into the variables passed by pointer
780 * into the constructor. If no neighbors are found, the output variables are
781 * not modified, i.e., the caller must initialize them.
787 * Initializes the action with given output locations.
789 * \param[out] closestPoint Index of the closest reference location.
790 * \param[out] minDist2 Minimum distance squared.
792 * The constructor call does not modify the pointed values, but only
793 * stores the pointers for later use.
794 * See the class description for additional semantics.
796 MindistAction(int *closestPoint, real *minDist2)
797 : closestPoint_(*closestPoint), minDist2_(*minDist2)
801 //! Processes a neighbor to find the nearest point.
802 bool operator()(int i, real r2)
816 GMX_DISALLOW_ASSIGN(MindistAction);
821 /********************************************************************
822 * AnalysisNeighborhood::Impl
825 class AnalysisNeighborhood::Impl
828 typedef AnalysisNeighborhoodSearch::ImplPointer SearchImplPointer;
829 typedef std::vector<SearchImplPointer> SearchList;
831 Impl() : cutoff_(0), excls_(NULL), mode_(eSearchMode_Automatic), bXY_(false)
836 SearchList::const_iterator i;
837 for (i = searchList_.begin(); i != searchList_.end(); ++i)
839 GMX_RELEASE_ASSERT(i->unique(),
840 "Dangling AnalysisNeighborhoodSearch reference");
844 SearchImplPointer getSearch();
846 tMPI::mutex createSearchMutex_;
847 SearchList searchList_;
849 const t_blocka *excls_;
854 AnalysisNeighborhood::Impl::SearchImplPointer
855 AnalysisNeighborhood::Impl::getSearch()
857 tMPI::lock_guard<tMPI::mutex> lock(createSearchMutex_);
858 // TODO: Consider whether this needs to/can be faster, e.g., by keeping a
859 // separate pool of unused search objects.
860 SearchList::const_iterator i;
861 for (i = searchList_.begin(); i != searchList_.end(); ++i)
868 SearchImplPointer search(new internal::AnalysisNeighborhoodSearchImpl(cutoff_));
869 searchList_.push_back(search);
873 /********************************************************************
874 * AnalysisNeighborhood
877 AnalysisNeighborhood::AnalysisNeighborhood()
882 AnalysisNeighborhood::~AnalysisNeighborhood()
886 void AnalysisNeighborhood::setCutoff(real cutoff)
888 GMX_RELEASE_ASSERT(impl_->searchList_.empty(),
889 "Changing the cutoff after initSearch() not currently supported");
890 impl_->cutoff_ = cutoff;
893 void AnalysisNeighborhood::setXYMode(bool bXY)
898 void AnalysisNeighborhood::setTopologyExclusions(const t_blocka *excls)
900 GMX_RELEASE_ASSERT(impl_->searchList_.empty(),
901 "Changing the exclusions after initSearch() not currently supported");
902 impl_->excls_ = excls;
905 void AnalysisNeighborhood::setMode(SearchMode mode)
910 AnalysisNeighborhood::SearchMode AnalysisNeighborhood::mode() const
915 AnalysisNeighborhoodSearch
916 AnalysisNeighborhood::initSearch(const t_pbc *pbc,
917 const AnalysisNeighborhoodPositions &positions)
919 Impl::SearchImplPointer search(impl_->getSearch());
920 search->init(mode(), impl_->bXY_, impl_->excls_, pbc, positions);
921 return AnalysisNeighborhoodSearch(search);
924 /********************************************************************
925 * AnalysisNeighborhoodSearch
928 AnalysisNeighborhoodSearch::AnalysisNeighborhoodSearch()
932 AnalysisNeighborhoodSearch::AnalysisNeighborhoodSearch(const ImplPointer &impl)
937 void AnalysisNeighborhoodSearch::reset()
942 AnalysisNeighborhood::SearchMode AnalysisNeighborhoodSearch::mode() const
944 GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
945 return (impl_->usesGridSearch()
946 ? AnalysisNeighborhood::eSearchMode_Grid
947 : AnalysisNeighborhood::eSearchMode_Simple);
950 bool AnalysisNeighborhoodSearch::isWithin(
951 const AnalysisNeighborhoodPositions &positions) const
953 GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
954 internal::AnalysisNeighborhoodPairSearchImpl pairSearch(*impl_);
955 pairSearch.startSearch(positions);
956 return pairSearch.searchNext(&withinAction);
959 real AnalysisNeighborhoodSearch::minimumDistance(
960 const AnalysisNeighborhoodPositions &positions) const
962 GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
963 internal::AnalysisNeighborhoodPairSearchImpl pairSearch(*impl_);
964 pairSearch.startSearch(positions);
965 real minDist2 = impl_->cutoffSquared();
966 int closestPoint = -1;
967 MindistAction action(&closestPoint, &minDist2);
968 (void)pairSearch.searchNext(action);
969 return sqrt(minDist2);
972 AnalysisNeighborhoodPair
973 AnalysisNeighborhoodSearch::nearestPoint(
974 const AnalysisNeighborhoodPositions &positions) const
976 GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
977 internal::AnalysisNeighborhoodPairSearchImpl pairSearch(*impl_);
978 pairSearch.startSearch(positions);
979 real minDist2 = impl_->cutoffSquared();
980 int closestPoint = -1;
981 MindistAction action(&closestPoint, &minDist2);
982 (void)pairSearch.searchNext(action);
983 return AnalysisNeighborhoodPair(closestPoint, 0, minDist2);
986 AnalysisNeighborhoodPairSearch
987 AnalysisNeighborhoodSearch::startPairSearch(
988 const AnalysisNeighborhoodPositions &positions) const
990 GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
991 Impl::PairSearchImplPointer pairSearch(impl_->getPairSearch());
992 pairSearch->startSearch(positions);
993 return AnalysisNeighborhoodPairSearch(pairSearch);
996 /********************************************************************
997 * AnalysisNeighborhoodPairSearch
1000 AnalysisNeighborhoodPairSearch::AnalysisNeighborhoodPairSearch(
1001 const ImplPointer &impl)
1006 bool AnalysisNeighborhoodPairSearch::findNextPair(AnalysisNeighborhoodPair *pair)
1008 bool bFound = impl_->searchNext(&withinAction);
1009 impl_->initFoundPair(pair);
1013 void AnalysisNeighborhoodPairSearch::skipRemainingPairsForTestPosition()
1015 impl_->nextTestPosition();