rvec_sub(maxBound, origin, size);
}
-} // namespace
+} // namespace
namespace internal
{
class AnalysisNeighborhoodSearchImpl
{
- public:
- typedef AnalysisNeighborhoodPairSearch::ImplPointer
- PairSearchImplPointer;
- typedef std::vector<PairSearchImplPointer> PairSearchList;
- typedef std::vector<std::vector<int> > CellList;
-
- explicit AnalysisNeighborhoodSearchImpl(real cutoff);
- ~AnalysisNeighborhoodSearchImpl();
-
- /*! \brief
- * Initializes the search with a given box and reference positions.
- *
- * \param[in] mode Search mode to use.
- * \param[in] bXY Whether to use 2D searching.
- * \param[in] excls Exclusions.
- * \param[in] pbc PBC information.
- * \param[in] positions Set of reference positions.
- */
- void init(AnalysisNeighborhood::SearchMode mode,
- bool bXY,
- const t_blocka *excls,
- const t_pbc *pbc,
- const AnalysisNeighborhoodPositions &positions);
- PairSearchImplPointer getPairSearch();
-
- real cutoffSquared() const { return cutoff2_; }
- bool usesGridSearch() const { return bGrid_; }
-
- private:
- /*! \brief
- * Determines a suitable grid size and sets up the cells.
- *
- * \param[in] box Box vectors (should not have zero vectors).
- * \param[in] bSingleCell If `true`, the corresponding dimension will
- * be forced to use a single cell.
- * \param[in] posCount Number of positions that will be put on the
- * grid.
- * \returns `false` if grid search is not suitable.
- */
- bool initGridCells(const matrix box, bool bSingleCell[DIM],
- int posCount);
- /*! \brief
- * Sets ua a search grid for a given box.
- *
- * \param[in] pbc Information about the box.
- * \param[in] posCount Number of positions in \p x.
- * \param[in] x Reference positions that will be put on the grid.
- * \param[in] bForce If `true`, grid searching will be used if at all
- * possible, even if a simple search might give better performance.
- * \returns `false` if grid search is not suitable.
- */
- bool initGrid(const t_pbc &pbc, int posCount, const rvec x[], bool bForce);
- /*! \brief
- * Maps a point into a grid cell.
- *
- * \param[in] x Point to map.
- * \param[out] cell Fractional cell coordinates of \p x on the grid.
- * \param[out] xout Coordinates to use.
- *
- * \p xout will be within the rectangular unit cell in dimensions where
- * the grid is periodic. For other dimensions, both \p xout and
- * \p cell can be outside the grid/unit cell.
- */
- void mapPointToGridCell(const rvec x, rvec cell, rvec xout) const;
- /*! \brief
- * Calculates linear index of a grid cell.
- *
- * \param[in] cell Cell indices (must be within the grid).
- * \returns Linear index of \p cell.
- */
- int getGridCellIndex(const ivec cell) const;
- /*! \brief
- * Calculates linear index of a grid cell from fractional coordinates.
- *
- * \param[in] cell Cell indices (must be within the grid).
- * \returns Linear index of \p cell.
- */
- int getGridCellIndex(const rvec cell) const;
- /*! \brief
- * Adds an index into a grid cell.
- *
- * \param[in] cell Fractional cell coordinates into which \p i should
- * be added.
- * \param[in] i Index to add.
- *
- * \p cell should satisfy the conditions that \p mapPointToGridCell()
- * produces.
- */
- void addToGridCell(const rvec cell, int i);
- /*! \brief
- * Initializes a cell pair loop for a dimension.
- *
- * \param[in] centerCell Fractional cell coordiates of the particle
- * for which pairs are being searched.
- * \param[in,out] cell Current/initial cell to loop over.
- * \param[in,out] upperBound Last cell to loop over.
- * \param[in] dim Dimension to initialize in this call.
- *
- * Initializes `cell[dim]` and `upperBound[dim]` for looping over
- * neighbors of a particle at position given by \p centerCell.
- * If 'dim != ZZ`, `cell[d]` (`d > dim`) set the plane/row of cells
- * for which the loop is initialized. The loop should then go from
- * `cell[dim]` until `upperBound[dim]`, inclusive.
- * `cell[d]` with `d < dim` or `upperBound[d]` with `d != dim` are not
- * modified by this function.
- *
- * `cell` and `upperBound` may be outside the grid for periodic
- * dimensions and need to be shifted separately: to simplify the
- * looping, the range is always (roughly) symmetric around the value in
- * `centerCell`.
- */
- void initCellRange(const rvec centerCell, ivec cell,
- ivec upperBound, int dim) const;
- /*! \brief
- * Computes the extent of the cutoff sphere on a particular cell edge.
- *
- * \param[in] centerCell Fractional cell coordiates of the particle
- * for which pairs are being searched.
- * \param[in] cell Current cell (for dimensions `>dim`).
- * \param[in] dim Dimension to compute in this call.
- * \returns Fractional extent of the cutoff sphere when looping
- * over cells in dimension `dim`, for `cell[d]` (`d > dim`).
- *
- * Input parameters are as for initCellRange(), except that if `cell`
- * is over a periodic boundary from `centerCell`, triclinic shifts
- * should have been applied to `centerCell` X/Y components.
- */
- real computeCutoffExtent(RVec centerCell, const ivec cell, int dim) const;
- /*! \brief
- * Advances cell pair loop to the next cell.
- *
- * \param[in] centerCell Fractional cell coordiates of the particle
- * for which pairs are being searched.
- * \param[in,out] cell Current (in)/next (out) cell in the loop.
- * \param[in,out] upperBound Last cell in the loop for each dimension.
- */
- bool nextCell(const rvec centerCell, ivec cell, ivec upperBound) const;
- /*! \brief
- * Calculates the index and shift of a grid cell during looping.
- *
- * \param[in] cell Unshifted cell index.
- * \param[out] shift Shift to apply to get the periodic distance
- * for distances between the cells.
- * \returns Grid cell index corresponding to `cell`.
- */
- int shiftCell(const ivec cell, rvec shift) const;
-
- //! Whether to try grid searching.
- bool bTryGrid_;
- //! The cutoff.
- real cutoff_;
- //! The cutoff squared.
- real cutoff2_;
- //! Whether to do searching in XY plane only.
- bool bXY_;
-
- //! Number of reference points for the current frame.
- int nref_;
- //! Reference point positions.
- const rvec *xref_;
- //! Reference position exclusion IDs.
- const int *refExclusionIds_;
- //! Reference position indices (NULL if no indices).
- const int *refIndices_;
- //! Exclusions.
- const t_blocka *excls_;
- //! PBC data.
- t_pbc pbc_;
-
- //! Whether grid searching is actually used for the current positions.
- bool bGrid_;
- //! false if the box is rectangular.
- bool bTric_;
- //! Whether the grid is periodic in a dimension.
- bool bGridPBC_[DIM];
- //! Array for storing in-unit-cell reference positions.
- std::vector<RVec> xrefAlloc_;
- //! Origin of the grid (zero for periodic dimensions).
- rvec gridOrigin_;
- //! Size of a single grid cell.
- rvec cellSize_;
- //! Inverse of \p cellSize_. Zero for dimensions where grid is not used.
- rvec invCellSize_;
- /*! \brief
- * Shift in cell coordinates (for triclinic boxes) in X when crossing
- * the Z periodic boundary.
- */
- real cellShiftZX_;
- /*! \brief
- * Shift in cell coordinates (for triclinic boxes) in Y when crossing
- * the Z periodic boundary.
- */
- real cellShiftZY_;
- /*! \brief
- * Shift in cell coordinates (for triclinic boxes) in X when crossing
- * the Y periodic boundary.
- */
- real cellShiftYX_;
- //! Number of cells along each dimension.
- ivec ncelldim_;
- //! Data structure to hold the grid cell contents.
- CellList cells_;
-
- Mutex createPairSearchMutex_;
- PairSearchList pairSearchList_;
-
- friend class AnalysisNeighborhoodPairSearchImpl;
-
- GMX_DISALLOW_COPY_AND_ASSIGN(AnalysisNeighborhoodSearchImpl);
+public:
+ typedef AnalysisNeighborhoodPairSearch::ImplPointer PairSearchImplPointer;
+ typedef std::vector<PairSearchImplPointer> PairSearchList;
+ typedef std::vector<std::vector<int>> CellList;
+
+ explicit AnalysisNeighborhoodSearchImpl(real cutoff);
+ ~AnalysisNeighborhoodSearchImpl();
+
+ /*! \brief
+ * Initializes the search with a given box and reference positions.
+ *
+ * \param[in] mode Search mode to use.
+ * \param[in] bXY Whether to use 2D searching.
+ * \param[in] excls Exclusions.
+ * \param[in] pbc PBC information.
+ * \param[in] positions Set of reference positions.
+ */
+ void init(AnalysisNeighborhood::SearchMode mode,
+ bool bXY,
+ const t_blocka* excls,
+ const t_pbc* pbc,
+ const AnalysisNeighborhoodPositions& positions);
+ PairSearchImplPointer getPairSearch();
+
+ real cutoffSquared() const { return cutoff2_; }
+ bool usesGridSearch() const { return bGrid_; }
+
+private:
+ /*! \brief
+ * Determines a suitable grid size and sets up the cells.
+ *
+ * \param[in] box Box vectors (should not have zero vectors).
+ * \param[in] bSingleCell If `true`, the corresponding dimension will
+ * be forced to use a single cell.
+ * \param[in] posCount Number of positions that will be put on the
+ * grid.
+ * \returns `false` if grid search is not suitable.
+ */
+ bool initGridCells(const matrix box, bool bSingleCell[DIM], int posCount);
+ /*! \brief
+ * Sets ua a search grid for a given box.
+ *
+ * \param[in] pbc Information about the box.
+ * \param[in] posCount Number of positions in \p x.
+ * \param[in] x Reference positions that will be put on the grid.
+ * \param[in] bForce If `true`, grid searching will be used if at all
+ * possible, even if a simple search might give better performance.
+ * \returns `false` if grid search is not suitable.
+ */
+ bool initGrid(const t_pbc& pbc, int posCount, const rvec x[], bool bForce);
+ /*! \brief
+ * Maps a point into a grid cell.
+ *
+ * \param[in] x Point to map.
+ * \param[out] cell Fractional cell coordinates of \p x on the grid.
+ * \param[out] xout Coordinates to use.
+ *
+ * \p xout will be within the rectangular unit cell in dimensions where
+ * the grid is periodic. For other dimensions, both \p xout and
+ * \p cell can be outside the grid/unit cell.
+ */
+ void mapPointToGridCell(const rvec x, rvec cell, rvec xout) const;
+ /*! \brief
+ * Calculates linear index of a grid cell.
+ *
+ * \param[in] cell Cell indices (must be within the grid).
+ * \returns Linear index of \p cell.
+ */
+ int getGridCellIndex(const ivec cell) const;
+ /*! \brief
+ * Calculates linear index of a grid cell from fractional coordinates.
+ *
+ * \param[in] cell Cell indices (must be within the grid).
+ * \returns Linear index of \p cell.
+ */
+ int getGridCellIndex(const rvec cell) const;
+ /*! \brief
+ * Adds an index into a grid cell.
+ *
+ * \param[in] cell Fractional cell coordinates into which \p i should
+ * be added.
+ * \param[in] i Index to add.
+ *
+ * \p cell should satisfy the conditions that \p mapPointToGridCell()
+ * produces.
+ */
+ void addToGridCell(const rvec cell, int i);
+ /*! \brief
+ * Initializes a cell pair loop for a dimension.
+ *
+ * \param[in] centerCell Fractional cell coordiates of the particle
+ * for which pairs are being searched.
+ * \param[in,out] cell Current/initial cell to loop over.
+ * \param[in,out] upperBound Last cell to loop over.
+ * \param[in] dim Dimension to initialize in this call.
+ *
+ * Initializes `cell[dim]` and `upperBound[dim]` for looping over
+ * neighbors of a particle at position given by \p centerCell.
+ * If 'dim != ZZ`, `cell[d]` (`d > dim`) set the plane/row of cells
+ * for which the loop is initialized. The loop should then go from
+ * `cell[dim]` until `upperBound[dim]`, inclusive.
+ * `cell[d]` with `d < dim` or `upperBound[d]` with `d != dim` are not
+ * modified by this function.
+ *
+ * `cell` and `upperBound` may be outside the grid for periodic
+ * dimensions and need to be shifted separately: to simplify the
+ * looping, the range is always (roughly) symmetric around the value in
+ * `centerCell`.
+ */
+ void initCellRange(const rvec centerCell, ivec cell, ivec upperBound, int dim) const;
+ /*! \brief
+ * Computes the extent of the cutoff sphere on a particular cell edge.
+ *
+ * \param[in] centerCell Fractional cell coordiates of the particle
+ * for which pairs are being searched.
+ * \param[in] cell Current cell (for dimensions `>dim`).
+ * \param[in] dim Dimension to compute in this call.
+ * \returns Fractional extent of the cutoff sphere when looping
+ * over cells in dimension `dim`, for `cell[d]` (`d > dim`).
+ *
+ * Input parameters are as for initCellRange(), except that if `cell`
+ * is over a periodic boundary from `centerCell`, triclinic shifts
+ * should have been applied to `centerCell` X/Y components.
+ */
+ real computeCutoffExtent(RVec centerCell, const ivec cell, int dim) const;
+ /*! \brief
+ * Advances cell pair loop to the next cell.
+ *
+ * \param[in] centerCell Fractional cell coordiates of the particle
+ * for which pairs are being searched.
+ * \param[in,out] cell Current (in)/next (out) cell in the loop.
+ * \param[in,out] upperBound Last cell in the loop for each dimension.
+ */
+ bool nextCell(const rvec centerCell, ivec cell, ivec upperBound) const;
+ /*! \brief
+ * Calculates the index and shift of a grid cell during looping.
+ *
+ * \param[in] cell Unshifted cell index.
+ * \param[out] shift Shift to apply to get the periodic distance
+ * for distances between the cells.
+ * \returns Grid cell index corresponding to `cell`.
+ */
+ int shiftCell(const ivec cell, rvec shift) const;
+
+ //! Whether to try grid searching.
+ bool bTryGrid_;
+ //! The cutoff.
+ real cutoff_;
+ //! The cutoff squared.
+ real cutoff2_;
+ //! Whether to do searching in XY plane only.
+ bool bXY_;
+
+ //! Number of reference points for the current frame.
+ int nref_;
+ //! Reference point positions.
+ const rvec* xref_;
+ //! Reference position exclusion IDs.
+ const int* refExclusionIds_;
+ //! Reference position indices (NULL if no indices).
+ const int* refIndices_;
+ //! Exclusions.
+ const t_blocka* excls_;
+ //! PBC data.
+ t_pbc pbc_;
+
+ //! Whether grid searching is actually used for the current positions.
+ bool bGrid_;
+ //! false if the box is rectangular.
+ bool bTric_;
+ //! Whether the grid is periodic in a dimension.
+ bool bGridPBC_[DIM];
+ //! Array for storing in-unit-cell reference positions.
+ std::vector<RVec> xrefAlloc_;
+ //! Origin of the grid (zero for periodic dimensions).
+ rvec gridOrigin_;
+ //! Size of a single grid cell.
+ rvec cellSize_;
+ //! Inverse of \p cellSize_. Zero for dimensions where grid is not used.
+ rvec invCellSize_;
+ /*! \brief
+ * Shift in cell coordinates (for triclinic boxes) in X when crossing
+ * the Z periodic boundary.
+ */
+ real cellShiftZX_;
+ /*! \brief
+ * Shift in cell coordinates (for triclinic boxes) in Y when crossing
+ * the Z periodic boundary.
+ */
+ real cellShiftZY_;
+ /*! \brief
+ * Shift in cell coordinates (for triclinic boxes) in X when crossing
+ * the Y periodic boundary.
+ */
+ real cellShiftYX_;
+ //! Number of cells along each dimension.
+ ivec ncelldim_;
+ //! Data structure to hold the grid cell contents.
+ CellList cells_;
+
+ Mutex createPairSearchMutex_;
+ PairSearchList pairSearchList_;
+
+ friend class AnalysisNeighborhoodPairSearchImpl;
+
+ GMX_DISALLOW_COPY_AND_ASSIGN(AnalysisNeighborhoodSearchImpl);
};
class AnalysisNeighborhoodPairSearchImpl
{
- public:
- explicit AnalysisNeighborhoodPairSearchImpl(const AnalysisNeighborhoodSearchImpl &search)
- : search_(search)
- {
- selfSearchMode_ = false;
- testPosCount_ = 0;
- testPositions_ = nullptr;
- testExclusionIds_ = nullptr;
- testIndices_ = nullptr;
- nexcl_ = 0;
- excl_ = nullptr;
- clear_rvec(xtest_);
- clear_rvec(testcell_);
- clear_ivec(currCell_);
- clear_ivec(cellBound_);
- reset(-1);
- }
+public:
+ explicit AnalysisNeighborhoodPairSearchImpl(const AnalysisNeighborhoodSearchImpl& search) :
+ search_(search)
+ {
+ selfSearchMode_ = false;
+ testPosCount_ = 0;
+ testPositions_ = nullptr;
+ testExclusionIds_ = nullptr;
+ testIndices_ = nullptr;
+ nexcl_ = 0;
+ excl_ = nullptr;
+ clear_rvec(xtest_);
+ clear_rvec(testcell_);
+ clear_ivec(currCell_);
+ clear_ivec(cellBound_);
+ reset(-1);
+ }
- //! Initializes a search to find reference positions neighboring \p x.
- void startSearch(const AnalysisNeighborhoodPositions &positions);
- //! Initializes a search to find reference position pairs.
- void startSelfSearch();
- //! Searches for the next neighbor.
- template <class Action>
- bool searchNext(Action action);
- //! Initializes a pair representing the pair found by searchNext().
- void initFoundPair(AnalysisNeighborhoodPair *pair) const;
- //! Advances to the next test position, skipping any remaining pairs.
- void nextTestPosition();
-
- private:
- //! Clears the loop indices.
- void reset(int testIndex);
- //! Checks whether a reference positiong should be excluded.
- bool isExcluded(int j);
-
- //! Parent search object.
- const AnalysisNeighborhoodSearchImpl &search_;
- //! Whether we are searching for ref-ref pairs.
- bool selfSearchMode_;
- //! Number of test positions.
- int testPosCount_;
- //! Reference to the test positions.
- const rvec *testPositions_;
- //! Reference to the test exclusion indices.
- const int *testExclusionIds_;
- //! Reference to the test position indices.
- const int *testIndices_;
- //! Number of excluded reference positions for current test particle.
- int nexcl_;
- //! Exclusions for current test particle.
- const int *excl_;
- //! Index of the currently active test position in \p testPositions_.
- int testIndex_;
- //! Stores test position during a pair loop.
- rvec xtest_;
- //! Stores the previous returned position during a pair loop.
- int previ_;
- //! Stores the pair distance corresponding to previ_.
- real prevr2_;
- //! Stores the shortest distance vector corresponding to previ_.
- rvec prevdx_;
- //! Stores the current exclusion index during loops.
- int exclind_;
- //! Stores the fractional test particle cell location during loops.
- rvec testcell_;
- //! Stores the cell index corresponding to testcell_.
- int testCellIndex_;
- //! Stores the current cell during pair loops.
- ivec currCell_;
- //! Stores the current loop upper bounds for each dimension during pair loops.
- ivec cellBound_;
- //! Stores the index within the current cell during pair loops.
- int prevcai_;
-
- GMX_DISALLOW_COPY_AND_ASSIGN(AnalysisNeighborhoodPairSearchImpl);
+ //! Initializes a search to find reference positions neighboring \p x.
+ void startSearch(const AnalysisNeighborhoodPositions& positions);
+ //! Initializes a search to find reference position pairs.
+ void startSelfSearch();
+ //! Searches for the next neighbor.
+ template<class Action>
+ bool searchNext(Action action);
+ //! Initializes a pair representing the pair found by searchNext().
+ void initFoundPair(AnalysisNeighborhoodPair* pair) const;
+ //! Advances to the next test position, skipping any remaining pairs.
+ void nextTestPosition();
+
+private:
+ //! Clears the loop indices.
+ void reset(int testIndex);
+ //! Checks whether a reference positiong should be excluded.
+ bool isExcluded(int j);
+
+ //! Parent search object.
+ const AnalysisNeighborhoodSearchImpl& search_;
+ //! Whether we are searching for ref-ref pairs.
+ bool selfSearchMode_;
+ //! Number of test positions.
+ int testPosCount_;
+ //! Reference to the test positions.
+ const rvec* testPositions_;
+ //! Reference to the test exclusion indices.
+ const int* testExclusionIds_;
+ //! Reference to the test position indices.
+ const int* testIndices_;
+ //! Number of excluded reference positions for current test particle.
+ int nexcl_;
+ //! Exclusions for current test particle.
+ const int* excl_;
+ //! Index of the currently active test position in \p testPositions_.
+ int testIndex_;
+ //! Stores test position during a pair loop.
+ rvec xtest_;
+ //! Stores the previous returned position during a pair loop.
+ int previ_;
+ //! Stores the pair distance corresponding to previ_.
+ real prevr2_;
+ //! Stores the shortest distance vector corresponding to previ_.
+ rvec prevdx_;
+ //! Stores the current exclusion index during loops.
+ int exclind_;
+ //! Stores the fractional test particle cell location during loops.
+ rvec testcell_;
+ //! Stores the cell index corresponding to testcell_.
+ int testCellIndex_;
+ //! Stores the current cell during pair loops.
+ ivec currCell_;
+ //! Stores the current loop upper bounds for each dimension during pair loops.
+ ivec cellBound_;
+ //! Stores the index within the current cell during pair loops.
+ int prevcai_;
+
+ GMX_DISALLOW_COPY_AND_ASSIGN(AnalysisNeighborhoodPairSearchImpl);
};
/********************************************************************
AnalysisNeighborhoodSearchImpl::AnalysisNeighborhoodSearchImpl(real cutoff)
{
- bTryGrid_ = true;
- cutoff_ = cutoff;
+ bTryGrid_ = true;
+ cutoff_ = cutoff;
if (cutoff_ <= 0)
{
- cutoff_ = cutoff2_ = GMX_REAL_MAX;
- bTryGrid_ = false;
+ cutoff_ = cutoff2_ = GMX_REAL_MAX;
+ bTryGrid_ = false;
}
else
{
- cutoff2_ = gmx::square(cutoff_);
+ cutoff2_ = gmx::square(cutoff_);
}
bXY_ = false;
nref_ = 0;
refIndices_ = nullptr;
std::memset(&pbc_, 0, sizeof(pbc_));
- bGrid_ = false;
- bTric_ = false;
- bGridPBC_[XX] = true;
- bGridPBC_[YY] = true;
- bGridPBC_[ZZ] = true;
+ bGrid_ = false;
+ bTric_ = false;
+ bGridPBC_[XX] = true;
+ bGridPBC_[YY] = true;
+ bGridPBC_[ZZ] = true;
clear_rvec(gridOrigin_);
clear_rvec(cellSize_);
PairSearchList::const_iterator i;
for (i = pairSearchList_.begin(); i != pairSearchList_.end(); ++i)
{
- GMX_RELEASE_ASSERT(i->unique(),
- "Dangling AnalysisNeighborhoodPairSearch reference");
+ GMX_RELEASE_ASSERT(i->unique(), "Dangling AnalysisNeighborhoodPairSearch reference");
}
}
-AnalysisNeighborhoodSearchImpl::PairSearchImplPointer
-AnalysisNeighborhoodSearchImpl::getPairSearch()
+AnalysisNeighborhoodSearchImpl::PairSearchImplPointer AnalysisNeighborhoodSearchImpl::getPairSearch()
{
lock_guard<Mutex> lock(createPairSearchMutex_);
// TODO: Consider whether this needs to/can be faster, e.g., by keeping a
return pairSearch;
}
-bool AnalysisNeighborhoodSearchImpl::initGridCells(
- const matrix box, bool bSingleCell[DIM], int posCount)
+bool AnalysisNeighborhoodSearchImpl::initGridCells(const matrix box, bool bSingleCell[DIM], int posCount)
{
// Determine the size of cubes where there are on average 10 positions.
// The loop takes care of cases where some of the box edges are shorter
{
break;
}
- targetsize = pow(volume * 10 / posCount, static_cast<real>(1./dimCount));
+ targetsize = pow(volume * 10 / posCount, static_cast<real>(1. / dimCount));
prevDimCount = dimCount;
}
}
}
totalCellCount *= cellCount;
- ncelldim_[dd] = cellCount;
+ ncelldim_[dd] = cellCount;
}
if (totalCellCount <= 3)
{
return true;
}
-bool AnalysisNeighborhoodSearchImpl::initGrid(
- const t_pbc &pbc, int posCount, const rvec x[], bool bForce)
+bool AnalysisNeighborhoodSearchImpl::initGrid(const t_pbc& pbc, int posCount, const rvec x[], bool bForce)
{
if (posCount == 0)
{
return false;
}
- bool bSingleCell[DIM] = {false, false, bXY_};
+ bool bSingleCell[DIM] = { false, false, bXY_ };
matrix box;
copy_mat(pbc.box, box);
// TODO: In principle, we could use the bounding box for periodic
// dimensions as well if the bounding box is sufficiently far from the box
// edges.
- rvec origin, boundingBoxSize;
+ rvec origin, boundingBoxSize;
computeBoundingBox(posCount, x, origin, boundingBoxSize);
clear_rvec(gridOrigin_);
for (int dd = 0; dd < DIM; ++dd)
// TODO: It could be better to avoid this when determining the cell
// size, but this can still remain here as a fallback to avoid
// incorrect results.
- if (std::ceil(2*cutoff_*invCellSize_[dd]) >= ncelldim_[dd])
+ if (std::ceil(2 * cutoff_ * invCellSize_[dd]) >= ncelldim_[dd])
{
// Cutoff is too close to half the box size for grid searching
// (it is not possible to find a single shift for every pair of
return true;
}
-void AnalysisNeighborhoodSearchImpl::mapPointToGridCell(const rvec x,
- rvec cell,
- rvec xout) const
+void AnalysisNeighborhoodSearchImpl::mapPointToGridCell(const rvec x, rvec cell, rvec xout) const
{
rvec xtmp;
rvec_sub(x, gridOrigin_, xtmp);
int AnalysisNeighborhoodSearchImpl::getGridCellIndex(const ivec cell) const
{
- GMX_ASSERT(cell[XX] >= 0 && cell[XX] < ncelldim_[XX],
- "Grid cell X index out of range");
- GMX_ASSERT(cell[YY] >= 0 && cell[YY] < ncelldim_[YY],
- "Grid cell Y index out of range");
- GMX_ASSERT(cell[ZZ] >= 0 && cell[ZZ] < ncelldim_[ZZ],
- "Grid cell Z index out of range");
- return cell[XX] + cell[YY] * ncelldim_[XX]
- + cell[ZZ] * ncelldim_[XX] * ncelldim_[YY];
+ GMX_ASSERT(cell[XX] >= 0 && cell[XX] < ncelldim_[XX], "Grid cell X index out of range");
+ GMX_ASSERT(cell[YY] >= 0 && cell[YY] < ncelldim_[YY], "Grid cell Y index out of range");
+ GMX_ASSERT(cell[ZZ] >= 0 && cell[ZZ] < ncelldim_[ZZ], "Grid cell Z index out of range");
+ return cell[XX] + cell[YY] * ncelldim_[XX] + cell[ZZ] * ncelldim_[XX] * ncelldim_[YY];
}
int AnalysisNeighborhoodSearchImpl::getGridCellIndex(const rvec cell) const
cells_[ci].push_back(i);
}
-void AnalysisNeighborhoodSearchImpl::initCellRange(
- const rvec centerCell, ivec currCell, ivec upperBound, int dim) const
+void AnalysisNeighborhoodSearchImpl::initCellRange(const rvec centerCell, ivec currCell, ivec upperBound, int dim) const
{
- RVec shiftedCenter(centerCell);
+ RVec shiftedCenter(centerCell);
// Shift the center to the cell coordinates of currCell, so that
// computeCutoffExtent() can assume simple rectangular grid.
if (bTric_)
upperBound[dim] = static_cast<int>(floor(endOffset));
}
-real AnalysisNeighborhoodSearchImpl::computeCutoffExtent(
- const RVec centerCell, const ivec cell, int dim) const
+real AnalysisNeighborhoodSearchImpl::computeCutoffExtent(const RVec centerCell, const ivec cell, int dim) const
{
if (dim == ZZ)
{
{
continue;
}
- dist2 += dimDist*dimDist*cellSize_[d]*cellSize_[d];
+ dist2 += dimDist * dimDist * cellSize_[d] * cellSize_[d];
}
if (dist2 >= cutoff2_)
{
return std::sqrt(cutoff2_ - dist2);
}
-bool AnalysisNeighborhoodSearchImpl::nextCell(
- const rvec centerCell, ivec cell, ivec upperBound) const
+bool AnalysisNeighborhoodSearchImpl::nextCell(const rvec centerCell, ivec cell, ivec upperBound) const
{
int dim = 0;
while (dim < DIM)
{
-next:
+ next:
++cell[dim];
if (cell[dim] > upperBound[dim])
{
return getGridCellIndex(shiftedCell);
}
-void AnalysisNeighborhoodSearchImpl::init(
- AnalysisNeighborhood::SearchMode mode,
- bool bXY,
- const t_blocka *excls,
- const t_pbc *pbc,
- const AnalysisNeighborhoodPositions &positions)
+void AnalysisNeighborhoodSearchImpl::init(AnalysisNeighborhood::SearchMode mode,
+ bool bXY,
+ const t_blocka* excls,
+ const t_pbc* pbc,
+ const AnalysisNeighborhoodPositions& positions)
{
GMX_RELEASE_ASSERT(positions.index_ == -1,
"Individual indexed positions not supported as reference");
{
if (pbc->ePBC != epbcXY && pbc->ePBC != epbcXYZ)
{
- std::string message =
- formatString("Computations in the XY plane are not supported with PBC type '%s'",
- epbc_names[pbc->ePBC]);
+ std::string message = formatString(
+ "Computations in the XY plane are not supported with PBC type '%s'",
+ epbc_names[pbc->ePBC]);
GMX_THROW(NotImplementedError(message));
}
- if (pbc->ePBC == epbcXYZ &&
- (std::fabs(pbc->box[ZZ][XX]) > GMX_REAL_EPS*pbc->box[ZZ][ZZ] ||
- std::fabs(pbc->box[ZZ][YY]) > GMX_REAL_EPS*pbc->box[ZZ][ZZ]))
+ if (pbc->ePBC == epbcXYZ
+ && (std::fabs(pbc->box[ZZ][XX]) > GMX_REAL_EPS * pbc->box[ZZ][ZZ]
+ || std::fabs(pbc->box[ZZ][YY]) > GMX_REAL_EPS * pbc->box[ZZ][ZZ]))
{
- GMX_THROW(NotImplementedError("Computations in the XY plane are not supported when the last box vector is not parallel to the Z axis"));
+ GMX_THROW(
+ NotImplementedError("Computations in the XY plane are not supported when the "
+ "last box vector is not parallel to the Z axis"));
}
// Use a single grid cell in Z direction.
matrix box;
previ_ = -1;
prevr2_ = 0.0;
clear_rvec(prevdx_);
- exclind_ = 0;
- prevcai_ = -1;
+ exclind_ = 0;
+ prevcai_ = -1;
if (testIndex_ >= 0 && testIndex_ < testPosCount_)
{
- const int index =
- (testIndices_ != nullptr ? testIndices_[testIndex] : testIndex);
+ const int index = (testIndices_ != nullptr ? testIndices_[testIndex] : testIndex);
if (search_.bGrid_)
{
search_.mapPointToGridCell(testPositions_[index], testcell_, xtest_);
}
if (search_.excls_ != nullptr)
{
- const int exclIndex = testExclusionIds_[index];
+ const int exclIndex = testExclusionIds_[index];
if (exclIndex < search_.excls_->nr)
{
const int startIndex = search_.excls_->index[exclIndex];
- nexcl_ = search_.excls_->index[exclIndex + 1] - startIndex;
- excl_ = &search_.excls_->a[startIndex];
+ nexcl_ = search_.excls_->index[exclIndex + 1] - startIndex;
+ excl_ = &search_.excls_->a[startIndex];
}
else
{
{
if (exclind_ < nexcl_)
{
- const int index =
- (search_.refIndices_ != nullptr ? search_.refIndices_[j] : j);
+ const int index = (search_.refIndices_ != nullptr ? search_.refIndices_[j] : j);
const int refId = search_.refExclusionIds_[index];
while (exclind_ < nexcl_ && excl_[exclind_] < refId)
{
return false;
}
-void AnalysisNeighborhoodPairSearchImpl::startSearch(
- const AnalysisNeighborhoodPositions &positions)
+void AnalysisNeighborhoodPairSearchImpl::startSearch(const AnalysisNeighborhoodPositions& positions)
{
selfSearchMode_ = false;
testPosCount_ = positions.count_;
reset(0);
}
-template <class Action>
+template<class Action>
bool AnalysisNeighborhoodPairSearchImpl::searchNext(Action action)
{
while (testIndex_ < testPosCount_)
do
{
rvec shift;
- const int ci = search_.shiftCell(currCell_, shift);
+ const int ci = search_.shiftCell(currCell_, shift);
if (selfSearchMode_ && ci > testCellIndex_)
{
continue;
{
continue;
}
- rvec dx;
+ rvec dx;
rvec_sub(search_.xref_[i], xtest_, dx);
rvec_sub(dx, shift, dx);
- const real r2
- = search_.bXY_
- ? dx[XX]*dx[XX] + dx[YY]*dx[YY]
- : norm2(dx);
+ const real r2 = search_.bXY_ ? dx[XX] * dx[XX] + dx[YY] * dx[YY] : norm2(dx);
if (r2 <= search_.cutoff2_)
{
if (action(i, r2, dx))
}
exclind_ = 0;
cai = 0;
- }
- while (search_.nextCell(testcell_, currCell_, cellBound_));
+ } while (search_.nextCell(testcell_, currCell_, cellBound_));
}
else
{
{
rvec_sub(search_.xref_[i], xtest_, dx);
}
- const real r2
- = search_.bXY_
- ? dx[XX]*dx[XX] + dx[YY]*dx[YY]
- : norm2(dx);
+ const real r2 = search_.bXY_ ? dx[XX] * dx[XX] + dx[YY] * dx[YY] : norm2(dx);
if (r2 <= search_.cutoff2_)
{
if (action(i, r2, dx))
return false;
}
-void AnalysisNeighborhoodPairSearchImpl::initFoundPair(
- AnalysisNeighborhoodPair *pair) const
+void AnalysisNeighborhoodPairSearchImpl::initFoundPair(AnalysisNeighborhoodPair* pair) const
{
if (previ_ < 0)
{
}
}
-} // namespace internal
+} // namespace internal
namespace
{
*/
class MindistAction
{
- public:
- /*! \brief
- * Initializes the action with given output locations.
- *
- * \param[out] closestPoint Index of the closest reference location.
- * \param[out] minDist2 Minimum distance squared.
- * \param[out] dx Shortest distance vector.
- *
- * The constructor call does not modify the pointed values, but only
- * stores the pointers for later use.
- * See the class description for additional semantics.
- */
- MindistAction(int *closestPoint, real *minDist2, rvec *dx) // NOLINT(readability-non-const-parameter)
- : closestPoint_(*closestPoint), minDist2_(*minDist2), dx_(*dx)
- {
- }
- //! Copies the action.
- MindistAction(const MindistAction &) = default;
+public:
+ /*! \brief
+ * Initializes the action with given output locations.
+ *
+ * \param[out] closestPoint Index of the closest reference location.
+ * \param[out] minDist2 Minimum distance squared.
+ * \param[out] dx Shortest distance vector.
+ *
+ * The constructor call does not modify the pointed values, but only
+ * stores the pointers for later use.
+ * See the class description for additional semantics.
+ */
+ MindistAction(int* closestPoint, real* minDist2, rvec* dx) // NOLINT(readability-non-const-parameter)
+ :
+ closestPoint_(*closestPoint),
+ minDist2_(*minDist2),
+ dx_(*dx)
+ {
+ }
+ //! Copies the action.
+ MindistAction(const MindistAction&) = default;
- //! Processes a neighbor to find the nearest point.
- bool operator()(int i, real r2, const rvec dx)
+ //! Processes a neighbor to find the nearest point.
+ bool operator()(int i, real r2, const rvec dx)
+ {
+ if (r2 < minDist2_)
{
- if (r2 < minDist2_)
- {
- closestPoint_ = i;
- minDist2_ = r2;
- copy_rvec(dx, dx_);
- }
- return false;
+ closestPoint_ = i;
+ minDist2_ = r2;
+ copy_rvec(dx, dx_);
}
+ return false;
+ }
- private:
- int &closestPoint_;
- real &minDist2_;
- rvec &dx_;
+private:
+ int& closestPoint_;
+ real& minDist2_;
+ rvec& dx_;
- GMX_DISALLOW_ASSIGN(MindistAction);
+ GMX_DISALLOW_ASSIGN(MindistAction);
};
-} // namespace
+} // namespace
/********************************************************************
* AnalysisNeighborhood::Impl
class AnalysisNeighborhood::Impl
{
- public:
- typedef AnalysisNeighborhoodSearch::ImplPointer SearchImplPointer;
- typedef std::vector<SearchImplPointer> SearchList;
+public:
+ typedef AnalysisNeighborhoodSearch::ImplPointer SearchImplPointer;
+ typedef std::vector<SearchImplPointer> SearchList;
- Impl()
- : cutoff_(0), excls_(nullptr), mode_(eSearchMode_Automatic), bXY_(false)
- {
- }
- ~Impl()
+ Impl() : cutoff_(0), excls_(nullptr), mode_(eSearchMode_Automatic), bXY_(false) {}
+ ~Impl()
+ {
+ SearchList::const_iterator i;
+ for (i = searchList_.begin(); i != searchList_.end(); ++i)
{
- SearchList::const_iterator i;
- for (i = searchList_.begin(); i != searchList_.end(); ++i)
- {
- GMX_RELEASE_ASSERT(i->unique(),
- "Dangling AnalysisNeighborhoodSearch reference");
- }
+ GMX_RELEASE_ASSERT(i->unique(), "Dangling AnalysisNeighborhoodSearch reference");
}
+ }
- SearchImplPointer getSearch();
+ SearchImplPointer getSearch();
- Mutex createSearchMutex_;
- SearchList searchList_;
- real cutoff_;
- const t_blocka *excls_;
- SearchMode mode_;
- bool bXY_;
+ Mutex createSearchMutex_;
+ SearchList searchList_;
+ real cutoff_;
+ const t_blocka* excls_;
+ SearchMode mode_;
+ bool bXY_;
};
-AnalysisNeighborhood::Impl::SearchImplPointer
-AnalysisNeighborhood::Impl::getSearch()
+AnalysisNeighborhood::Impl::SearchImplPointer AnalysisNeighborhood::Impl::getSearch()
{
lock_guard<Mutex> lock(createSearchMutex_);
// TODO: Consider whether this needs to/can be faster, e.g., by keeping a
* AnalysisNeighborhood
*/
-AnalysisNeighborhood::AnalysisNeighborhood()
- : impl_(new Impl)
-{
-}
+AnalysisNeighborhood::AnalysisNeighborhood() : impl_(new Impl) {}
-AnalysisNeighborhood::~AnalysisNeighborhood()
-{
-}
+AnalysisNeighborhood::~AnalysisNeighborhood() {}
void AnalysisNeighborhood::setCutoff(real cutoff)
{
impl_->bXY_ = bXY;
}
-void AnalysisNeighborhood::setTopologyExclusions(const t_blocka *excls)
+void AnalysisNeighborhood::setTopologyExclusions(const t_blocka* excls)
{
GMX_RELEASE_ASSERT(impl_->searchList_.empty(),
"Changing the exclusions after initSearch() not currently supported");
return impl_->mode_;
}
-AnalysisNeighborhoodSearch
-AnalysisNeighborhood::initSearch(const t_pbc *pbc,
- const AnalysisNeighborhoodPositions &positions)
+AnalysisNeighborhoodSearch AnalysisNeighborhood::initSearch(const t_pbc* pbc,
+ const AnalysisNeighborhoodPositions& positions)
{
Impl::SearchImplPointer search(impl_->getSearch());
- search->init(mode(), impl_->bXY_, impl_->excls_,
- pbc, positions);
+ search->init(mode(), impl_->bXY_, impl_->excls_, pbc, positions);
return AnalysisNeighborhoodSearch(search);
}
* AnalysisNeighborhoodSearch
*/
-AnalysisNeighborhoodSearch::AnalysisNeighborhoodSearch()
-{
-}
+AnalysisNeighborhoodSearch::AnalysisNeighborhoodSearch() {}
-AnalysisNeighborhoodSearch::AnalysisNeighborhoodSearch(const ImplPointer &impl)
- : impl_(impl)
-{
-}
+AnalysisNeighborhoodSearch::AnalysisNeighborhoodSearch(const ImplPointer& impl) : impl_(impl) {}
void AnalysisNeighborhoodSearch::reset()
{
AnalysisNeighborhood::SearchMode AnalysisNeighborhoodSearch::mode() const
{
GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
- return (impl_->usesGridSearch()
- ? AnalysisNeighborhood::eSearchMode_Grid
- : AnalysisNeighborhood::eSearchMode_Simple);
+ return (impl_->usesGridSearch() ? AnalysisNeighborhood::eSearchMode_Grid
+ : AnalysisNeighborhood::eSearchMode_Simple);
}
-bool AnalysisNeighborhoodSearch::isWithin(
- const AnalysisNeighborhoodPositions &positions) const
+bool AnalysisNeighborhoodSearch::isWithin(const AnalysisNeighborhoodPositions& positions) const
{
GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
internal::AnalysisNeighborhoodPairSearchImpl pairSearch(*impl_);
return pairSearch.searchNext(&withinAction);
}
-real AnalysisNeighborhoodSearch::minimumDistance(
- const AnalysisNeighborhoodPositions &positions) const
+real AnalysisNeighborhoodSearch::minimumDistance(const AnalysisNeighborhoodPositions& positions) const
{
GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
internal::AnalysisNeighborhoodPairSearchImpl pairSearch(*impl_);
pairSearch.startSearch(positions);
real minDist2 = impl_->cutoffSquared();
int closestPoint = -1;
- rvec dx = {0.0, 0.0, 0.0};
+ rvec dx = { 0.0, 0.0, 0.0 };
MindistAction action(&closestPoint, &minDist2, &dx);
(void)pairSearch.searchNext(action);
return std::sqrt(minDist2);
}
-AnalysisNeighborhoodPair
-AnalysisNeighborhoodSearch::nearestPoint(
- const AnalysisNeighborhoodPositions &positions) const
+AnalysisNeighborhoodPair AnalysisNeighborhoodSearch::nearestPoint(const AnalysisNeighborhoodPositions& positions) const
{
GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
internal::AnalysisNeighborhoodPairSearchImpl pairSearch(*impl_);
pairSearch.startSearch(positions);
real minDist2 = impl_->cutoffSquared();
int closestPoint = -1;
- rvec dx = {0.0, 0.0, 0.0};
+ rvec dx = { 0.0, 0.0, 0.0 };
MindistAction action(&closestPoint, &minDist2, &dx);
(void)pairSearch.searchNext(action);
return AnalysisNeighborhoodPair(closestPoint, 0, minDist2, dx);
}
-AnalysisNeighborhoodPairSearch
-AnalysisNeighborhoodSearch::startSelfPairSearch() const
+AnalysisNeighborhoodPairSearch AnalysisNeighborhoodSearch::startSelfPairSearch() const
{
GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
Impl::PairSearchImplPointer pairSearch(impl_->getPairSearch());
}
AnalysisNeighborhoodPairSearch
-AnalysisNeighborhoodSearch::startPairSearch(
- const AnalysisNeighborhoodPositions &positions) const
+AnalysisNeighborhoodSearch::startPairSearch(const AnalysisNeighborhoodPositions& positions) const
{
GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
Impl::PairSearchImplPointer pairSearch(impl_->getPairSearch());
* AnalysisNeighborhoodPairSearch
*/
-AnalysisNeighborhoodPairSearch::AnalysisNeighborhoodPairSearch(
- const ImplPointer &impl)
- : impl_(impl)
+AnalysisNeighborhoodPairSearch::AnalysisNeighborhoodPairSearch(const ImplPointer& impl) :
+ impl_(impl)
{
}
-bool AnalysisNeighborhoodPairSearch::findNextPair(AnalysisNeighborhoodPair *pair)
+bool AnalysisNeighborhoodPairSearch::findNextPair(AnalysisNeighborhoodPair* pair)
{
bool bFound = impl_->searchNext(&withinAction);
impl_->initFoundPair(pair);
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff