src/gromacs/selection/nbsearch.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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.
   8  *
   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.
  13  *
  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.
  18  *
  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.
  23  *
  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.
  31  *
  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.
  34  */
  35 /*! \internal \file
  36  * \brief
  37  * Implements neighborhood searching for analysis (from nbsearch.h).
  38  *
  39  * \todo
  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
  42  *     handling easier.
  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.
  51  *
  52  * \author Teemu Murtola <teemu.murtola@gmail.com>
  53  * \ingroup module_selection
  54  */
  55 #include "gmxpre.h"
  56
  57 #include "nbsearch.h"
  58
  59 #include <cmath>
  60 #include <cstring>
  61
  62 #include <algorithm>
  63 #include <vector>
  64
  65 #include "thread_mpi/mutex.h"
  66
  67 #include "gromacs/legacyheaders/names.h"
  68 #include "gromacs/math/vec.h"
  69 #include "gromacs/pbcutil/pbc.h"
  70 #include "gromacs/selection/position.h"
  71 #include "gromacs/topology/block.h"
  72 #include "gromacs/utility/arrayref.h"
  73 #include "gromacs/utility/exceptions.h"
  74 #include "gromacs/utility/gmxassert.h"
  75 #include "gromacs/utility/smalloc.h"
  76 #include "gromacs/utility/stringutil.h"
  77
  78 namespace gmx
  79 {
  80
  81 namespace internal
  82 {
  83
  84 /********************************************************************
  85  * Implementation class declarations
  86  */
  87
  88 class AnalysisNeighborhoodSearchImpl
  89 {
  90     public:
  91         typedef AnalysisNeighborhoodPairSearch::ImplPointer
  92             PairSearchImplPointer;
  93         typedef std::vector<PairSearchImplPointer> PairSearchList;
  94         typedef std::vector<std::vector<int> > CellList;
  95
  96         explicit AnalysisNeighborhoodSearchImpl(real cutoff);
  97         ~AnalysisNeighborhoodSearchImpl();
  98
  99         /*! \brief
 100          * Initializes the search with a given box and reference positions.
 101          *
 102          * \param[in]     mode      Search mode to use.
 103          * \param[in]     bXY       Whether to use 2D searching.
 104          * \param[in]     excls     Exclusions.
 105          * \param[in]     pbc       PBC information.
 106          * \param[in]     positions Set of reference positions.
 107          */
 108         void init(AnalysisNeighborhood::SearchMode     mode,
 109                   bool                                 bXY,
 110                   const t_blocka                      *excls,
 111                   const t_pbc                         *pbc,
 112                   const AnalysisNeighborhoodPositions &positions);
 113         PairSearchImplPointer getPairSearch();
 114
 115         real cutoffSquared() const { return cutoff2_; }
 116         bool usesGridSearch() const { return bGrid_; }
 117
 118     private:
 119         //! Calculates offsets to neighboring grid cells that should be considered.
 120         void initGridCellNeighborList();
 121         /*! \brief
 122          * Determines a suitable grid size and sets up the cells.
 123          *
 124          * \param[in]     pbc  Information about the box.
 125          * \returns  false if grid search is not suitable.
 126          */
 127         bool initGridCells(const t_pbc &pbc);
 128         /*! \brief
 129          * Sets ua a search grid for a given box.
 130          *
 131          * \param[in]     pbc  Information about the box.
 132          * \returns  false if grid search is not suitable.
 133          */
 134         bool initGrid(const t_pbc &pbc);
 135         /*! \brief
 136          * Maps a point into a grid cell.
 137          *
 138          * \param[in]  x    Point to map.
 139          * \param[out] cell Indices of the grid cell in which \p x lies.
 140          * \param[out] xout Coordinates to use
 141          *     (will be within the triclinic unit cell).
 142          */
 143         void mapPointToGridCell(const rvec x, ivec cell, rvec xout) const;
 144         /*! \brief
 145          * Calculates linear index of a grid cell.
 146          *
 147          * \param[in]  cell Cell indices.
 148          * \returns    Linear index of \p cell.
 149          */
 150         int getGridCellIndex(const ivec cell) const;
 151         /*! \brief
 152          * Adds an index into a grid cell.
 153          *
 154          * \param[in]     cell Cell into which \p i should be added.
 155          * \param[in]     i    Index to add.
 156          */
 157         void addToGridCell(const ivec cell, int i);
 158         /*! \brief
 159          * Calculates the index of a neighboring grid cell.
 160          *
 161          * \param[in]  sourceCell Location of the source cell.
 162          * \param[in]  index      Index of the neighbor to calculate.
 163          * \param[out] shift      Shift to apply to get the periodic distance
 164          *     for distances between the cells.
 165          * \returns    Grid cell index of the neighboring cell.
 166          */
 167         int getNeighboringCell(const ivec sourceCell, int index, rvec shift) const;
 168
 169         //! Whether to try grid searching.
 170         bool                    bTryGrid_;
 171         //! The cutoff.
 172         real                    cutoff_;
 173         //! The cutoff squared.
 174         real                    cutoff2_;
 175         //! Whether to do searching in XY plane only.
 176         bool                    bXY_;
 177
 178         //! Number of reference points for the current frame.
 179         int                     nref_;
 180         //! Reference point positions.
 181         const rvec             *xref_;
 182         //! Reference position exclusion IDs.
 183         const int              *refExclusionIds_;
 184         //! Exclusions.
 185         const t_blocka         *excls_;
 186         //! PBC data.
 187         t_pbc                   pbc_;
 188
 189         //! Whether grid searching is actually used for the current positions.
 190         bool                    bGrid_;
 191         //! Array allocated for storing in-unit-cell reference positions.
 192         rvec                   *xref_alloc_;
 193         //! Allocation count for xref_alloc.
 194         int                     xref_nalloc_;
 195         //! false if the box is rectangular.
 196         bool                    bTric_;
 197         //! Box vectors of a single grid cell.
 198         matrix                  cellbox_;
 199         //! The reciprocal cell vectors as columns; the inverse of \p cellbox.
 200         matrix                  recipcell_;
 201         //! Number of cells along each dimension.
 202         ivec                    ncelldim_;
 203         //! Data structure to hold the grid cell contents.
 204         CellList                cells_;
 205         //! Number of neighboring cells to consider.
 206         int                     ngridnb_;
 207         //! Offsets of the neighboring cells to consider.
 208         ivec                   *gnboffs_;
 209         //! Allocation count for \p gnboffs.
 210         int                     gnboffs_nalloc_;
 211
 212         tMPI::mutex             createPairSearchMutex_;
 213         PairSearchList          pairSearchList_;
 214
 215         friend class AnalysisNeighborhoodPairSearchImpl;
 216
 217         GMX_DISALLOW_COPY_AND_ASSIGN(AnalysisNeighborhoodSearchImpl);
 218 };
 219
 220 class AnalysisNeighborhoodPairSearchImpl
 221 {
 222     public:
 223         explicit AnalysisNeighborhoodPairSearchImpl(const AnalysisNeighborhoodSearchImpl &search)
 224             : search_(search)
 225         {
 226             testExclusionIds_ = NULL;
 227             nexcl_            = 0;
 228             excl_             = NULL;
 229             clear_rvec(xtest_);
 230             clear_ivec(testcell_);
 231             reset(-1);
 232         }
 233
 234         //! Initializes a search to find reference positions neighboring \p x.
 235         void startSearch(const AnalysisNeighborhoodPositions &positions);
 236         //! Searches for the next neighbor.
 237         template <class Action>
 238         bool searchNext(Action action);
 239         //! Initializes a pair representing the pair found by searchNext().
 240         void initFoundPair(AnalysisNeighborhoodPair *pair) const;
 241         //! Advances to the next test position, skipping any remaining pairs.
 242         void nextTestPosition();
 243
 244     private:
 245         //! Clears the loop indices.
 246         void reset(int testIndex);
 247         //! Checks whether a reference positiong should be excluded.
 248         bool isExcluded(int j);
 249
 250         //! Parent search object.
 251         const AnalysisNeighborhoodSearchImpl   &search_;
 252         //! Reference to the test positions.
 253         ConstArrayRef<rvec>                     testPositions_;
 254         //! Reference to the test exclusion indices.
 255         const int                              *testExclusionIds_;
 256         //! Number of excluded reference positions for current test particle.
 257         int                                     nexcl_;
 258         //! Exclusions for current test particle.
 259         const int                              *excl_;
 260         //! Index of the currently active test position in \p testPositions_.
 261         int                                     testIndex_;
 262         //! Stores test position during a pair loop.
 263         rvec                                    xtest_;
 264         //! Stores the previous returned position during a pair loop.
 265         int                                     previ_;
 266         //! Stores the pair distance corresponding to previ_;
 267         real                                    prevr2_;
 268         //! Stores the current exclusion index during loops.
 269         int                                     exclind_;
 270         //! Stores the test particle cell index during loops.
 271         ivec                                    testcell_;
 272         //! Stores the current cell neighbor index during pair loops.
 273         int                                     prevnbi_;
 274         //! Stores the index within the current cell during pair loops.
 275         int                                     prevcai_;
 276
 277         GMX_DISALLOW_COPY_AND_ASSIGN(AnalysisNeighborhoodPairSearchImpl);
 278 };
 279
 280 /********************************************************************
 281  * AnalysisNeighborhoodSearchImpl
 282  */
 283
 284 AnalysisNeighborhoodSearchImpl::AnalysisNeighborhoodSearchImpl(real cutoff)
 285 {
 286     bTryGrid_       = true;
 287     cutoff_         = cutoff;
 288     if (cutoff_ <= 0)
 289     {
 290         cutoff_     = GMX_REAL_MAX;
 291         bTryGrid_   = false;
 292     }
 293     cutoff2_        = sqr(cutoff_);
 294     bXY_            = false;
 295
 296     nref_            = 0;
 297     xref_            = NULL;
 298     refExclusionIds_ = NULL;
 299     std::memset(&pbc_, 0, sizeof(pbc_));
 300
 301     bGrid_          = false;
 302
 303     xref_alloc_     = NULL;
 304     xref_nalloc_    = 0;
 305     bTric_          = false;
 306     clear_mat(cellbox_);
 307     clear_mat(recipcell_);
 308     clear_ivec(ncelldim_);
 309
 310     ngridnb_        = 0;
 311     gnboffs_        = NULL;
 312     gnboffs_nalloc_ = 0;
 313 }
 314
 315 AnalysisNeighborhoodSearchImpl::~AnalysisNeighborhoodSearchImpl()
 316 {
 317     PairSearchList::const_iterator i;
 318     for (i = pairSearchList_.begin(); i != pairSearchList_.end(); ++i)
 319     {
 320         GMX_RELEASE_ASSERT(i->unique(),
 321                            "Dangling AnalysisNeighborhoodPairSearch reference");
 322     }
 323     sfree(xref_alloc_);
 324     sfree(gnboffs_);
 325 }
 326
 327 AnalysisNeighborhoodSearchImpl::PairSearchImplPointer
 328 AnalysisNeighborhoodSearchImpl::getPairSearch()
 329 {
 330     tMPI::lock_guard<tMPI::mutex> lock(createPairSearchMutex_);
 331     // TODO: Consider whether this needs to/can be faster, e.g., by keeping a
 332     // separate pool of unused search objects.
 333     PairSearchList::const_iterator i;
 334     for (i = pairSearchList_.begin(); i != pairSearchList_.end(); ++i)
 335     {
 336         if (i->unique())
 337         {
 338             return *i;
 339         }
 340     }
 341     PairSearchImplPointer pairSearch(new AnalysisNeighborhoodPairSearchImpl(*this));
 342     pairSearchList_.push_back(pairSearch);
 343     return pairSearch;
 344 }
 345
 346 void AnalysisNeighborhoodSearchImpl::initGridCellNeighborList()
 347 {
 348     int   maxx, maxy, maxz;
 349     real  rvnorm;
 350
 351     /* Find the extent of the sphere in triclinic coordinates */
 352     maxz   = static_cast<int>(cutoff_ * recipcell_[ZZ][ZZ]) + 1;
 353     rvnorm = sqrt(sqr(recipcell_[YY][YY]) + sqr(recipcell_[ZZ][YY]));
 354     maxy   = static_cast<int>(cutoff_ * rvnorm) + 1;
 355     rvnorm = sqrt(sqr(recipcell_[XX][XX]) + sqr(recipcell_[YY][XX])
 356                   + sqr(recipcell_[ZZ][XX]));
 357     maxx   = static_cast<int>(cutoff_ * rvnorm) + 1;
 358
 359     /* Calculate the number of cells and reallocate if necessary */
 360     ngridnb_ = (2 * maxx + 1) * (2 * maxy + 1) * (2 * maxz + 1);
 361     if (gnboffs_nalloc_ < ngridnb_)
 362     {
 363         gnboffs_nalloc_ = ngridnb_;
 364         srenew(gnboffs_, gnboffs_nalloc_);
 365     }
 366
 367     /* Store the whole cube */
 368     /* TODO: Prune off corners that are not needed */
 369     int i = 0;
 370     for (int x = -maxx; x <= maxx; ++x)
 371     {
 372         for (int y = -maxy; y <= maxy; ++y)
 373         {
 374             for (int z = -maxz; z <= maxz; ++z)
 375             {
 376                 gnboffs_[i][XX] = x;
 377                 gnboffs_[i][YY] = y;
 378                 gnboffs_[i][ZZ] = z;
 379                 ++i;
 380             }
 381         }
 382     }
 383 }
 384
 385 bool AnalysisNeighborhoodSearchImpl::initGridCells(const t_pbc &pbc)
 386 {
 387     const real targetsize =
 388         pow(pbc.box[XX][XX] * pbc.box[YY][YY] * pbc.box[ZZ][ZZ]
 389             * 10 / nref_, static_cast<real>(1./3.));
 390
 391     int cellCount = 1;
 392     for (int dd = 0; dd < DIM; ++dd)
 393     {
 394         ncelldim_[dd] = static_cast<int>(pbc.box[dd][dd] / targetsize);
 395         cellCount    *= ncelldim_[dd];
 396         if (ncelldim_[dd] < 3)
 397         {
 398             return false;
 399         }
 400     }
 401     // Never decrease the size of the cell vector to avoid reallocating
 402     // memory for the nested vectors.  The actual size of the vector is not
 403     // used outside this function.
 404     if (cells_.size() < static_cast<size_t>(cellCount))
 405     {
 406         cells_.resize(cellCount);
 407     }
 408     for (int ci = 0; ci < cellCount; ++ci)
 409     {
 410         cells_[ci].clear();
 411     }
 412     return true;
 413 }
 414
 415 bool AnalysisNeighborhoodSearchImpl::initGrid(const t_pbc &pbc)
 416 {
 417     /* TODO: This check could be improved. */
 418     if (0.5*pbc.max_cutoff2 < cutoff2_)
 419     {
 420         return false;
 421     }
 422
 423     if (!initGridCells(pbc))
 424     {
 425         return false;
 426     }
 427
 428     bTric_ = TRICLINIC(pbc.box);
 429     if (bTric_)
 430     {
 431         for (int dd = 0; dd < DIM; ++dd)
 432         {
 433             svmul(1.0 / ncelldim_[dd], pbc.box[dd], cellbox_[dd]);
 434         }
 435         m_inv_ur0(cellbox_, recipcell_);
 436     }
 437     else
 438     {
 439         for (int dd = 0; dd < DIM; ++dd)
 440         {
 441             cellbox_[dd][dd]   = pbc.box[dd][dd] / ncelldim_[dd];
 442             recipcell_[dd][dd] = 1.0 / cellbox_[dd][dd];
 443         }
 444     }
 445     initGridCellNeighborList();
 446     return true;
 447 }
 448
 449 void AnalysisNeighborhoodSearchImpl::mapPointToGridCell(const rvec x,
 450                                                         ivec       cell,
 451                                                         rvec       xout) const
 452 {
 453     rvec xtmp;
 454     copy_rvec(x, xtmp);
 455     if (bTric_)
 456     {
 457         rvec tx;
 458         tmvmul_ur0(recipcell_, xtmp, tx);
 459         for (int dd = 0; dd < DIM; ++dd)
 460         {
 461             const int cellCount = ncelldim_[dd];
 462             int       cellIndex = static_cast<int>(floor(tx[dd]));
 463             while (cellIndex < 0)
 464             {
 465                 cellIndex += cellCount;
 466                 rvec_add(xtmp, pbc_.box[dd], xtmp);
 467             }
 468             while (cellIndex >= cellCount)
 469             {
 470                 cellIndex -= cellCount;
 471                 rvec_sub(xtmp, pbc_.box[dd], xtmp);
 472             }
 473             cell[dd] = cellIndex;
 474         }
 475     }
 476     else
 477     {
 478         for (int dd = 0; dd < DIM; ++dd)
 479         {
 480             const int cellCount = ncelldim_[dd];
 481             int       cellIndex = static_cast<int>(floor(xtmp[dd] * recipcell_[dd][dd]));
 482             while (cellIndex < 0)
 483             {
 484                 cellIndex += cellCount;
 485                 xtmp[dd]  += pbc_.box[dd][dd];
 486             }
 487             while (cellIndex >= cellCount)
 488             {
 489                 cellIndex -= cellCount;
 490                 xtmp[dd]  -= pbc_.box[dd][dd];
 491             }
 492             cell[dd] = cellIndex;
 493         }
 494     }
 495     copy_rvec(xtmp, xout);
 496 }
 497
 498 int AnalysisNeighborhoodSearchImpl::getGridCellIndex(const ivec cell) const
 499 {
 500     GMX_ASSERT(cell[XX] >= 0 && cell[XX] < ncelldim_[XX],
 501                "Grid cell X index out of range");
 502     GMX_ASSERT(cell[YY] >= 0 && cell[YY] < ncelldim_[YY],
 503                "Grid cell Y index out of range");
 504     GMX_ASSERT(cell[ZZ] >= 0 && cell[ZZ] < ncelldim_[ZZ],
 505                "Grid cell Z index out of range");
 506     return cell[XX] + cell[YY] * ncelldim_[XX]
 507            + cell[ZZ] * ncelldim_[XX] * ncelldim_[YY];
 508 }
 509
 510 void AnalysisNeighborhoodSearchImpl::addToGridCell(const ivec cell, int i)
 511 {
 512     const int ci = getGridCellIndex(cell);
 513     cells_[ci].push_back(i);
 514 }
 515
 516 int AnalysisNeighborhoodSearchImpl::getNeighboringCell(
 517         const ivec sourceCell, int index, rvec shift) const
 518 {
 519     ivec cell;
 520     ivec_add(sourceCell, gnboffs_[index], cell);
 521
 522     // TODO: Consider unifying with the similar shifting code in
 523     // mapPointToGridCell().
 524     clear_rvec(shift);
 525     for (int d = 0; d < DIM; ++d)
 526     {
 527         const int cellCount = ncelldim_[d];
 528         if (cell[d] < 0)
 529         {
 530             cell[d] += cellCount;
 531             rvec_add(shift, pbc_.box[d], shift);
 532         }
 533         else if (cell[d] >= cellCount)
 534         {
 535             cell[d] -= cellCount;
 536             rvec_sub(shift, pbc_.box[d], shift);
 537         }
 538     }
 539
 540     return getGridCellIndex(cell);
 541 }
 542
 543 void AnalysisNeighborhoodSearchImpl::init(
 544         AnalysisNeighborhood::SearchMode     mode,
 545         bool                                 bXY,
 546         const t_blocka                      *excls,
 547         const t_pbc                         *pbc,
 548         const AnalysisNeighborhoodPositions &positions)
 549 {
 550     GMX_RELEASE_ASSERT(positions.index_ == -1,
 551                        "Individual indexed positions not supported as reference");
 552     bXY_ = bXY;
 553     if (bXY_ && pbc->ePBC != epbcNONE)
 554     {
 555         if (pbc->ePBC != epbcXY && pbc->ePBC != epbcXYZ)
 556         {
 557             std::string message =
 558                 formatString("Computations in the XY plane are not supported with PBC type '%s'",
 559                              EPBC(pbc->ePBC));
 560             GMX_THROW(NotImplementedError(message));
 561         }
 562         if (std::fabs(pbc->box[ZZ][XX]) > GMX_REAL_EPS*pbc->box[ZZ][ZZ] ||
 563             std::fabs(pbc->box[ZZ][YY]) > GMX_REAL_EPS*pbc->box[ZZ][ZZ])
 564         {
 565             GMX_THROW(NotImplementedError("Computations in the XY plane are not supported when the last box vector is not parallel to the Z axis"));
 566         }
 567         set_pbc(&pbc_, epbcXY, const_cast<rvec *>(pbc->box));
 568     }
 569     else if (pbc != NULL)
 570     {
 571         pbc_  = *pbc;
 572     }
 573     else
 574     {
 575         pbc_.ePBC = epbcNONE;
 576     }
 577     nref_ = positions.count_;
 578     // TODO: Consider whether it would be possible to support grid searching in
 579     // more cases.
 580     if (mode == AnalysisNeighborhood::eSearchMode_Simple
 581         || pbc_.ePBC != epbcXYZ)
 582     {
 583         bGrid_ = false;
 584     }
 585     else if (bTryGrid_)
 586     {
 587         // TODO: Actually implement forcing eSearchMode_Grid
 588         bGrid_ = initGrid(pbc_);
 589     }
 590     if (bGrid_)
 591     {
 592         if (xref_nalloc_ < nref_)
 593         {
 594             srenew(xref_alloc_, nref_);
 595             xref_nalloc_ = nref_;
 596         }
 597         xref_ = xref_alloc_;
 598
 599         for (int i = 0; i < nref_; ++i)
 600         {
 601             ivec refcell;
 602             mapPointToGridCell(positions.x_[i], refcell, xref_alloc_[i]);
 603             addToGridCell(refcell, i);
 604         }
 605     }
 606     else
 607     {
 608         xref_ = positions.x_;
 609     }
 610     excls_           = excls;
 611     refExclusionIds_ = NULL;
 612     if (excls != NULL)
 613     {
 614         // TODO: Check that the IDs are ascending, or remove the limitation.
 615         refExclusionIds_ = positions.exclusionIds_;
 616         GMX_RELEASE_ASSERT(refExclusionIds_ != NULL,
 617                            "Exclusion IDs must be set for reference positions "
 618                            "when exclusions are enabled");
 619     }
 620 }
 621
 622 /********************************************************************
 623  * AnalysisNeighborhoodPairSearchImpl
 624  */
 625
 626 void AnalysisNeighborhoodPairSearchImpl::reset(int testIndex)
 627 {
 628     testIndex_ = testIndex;
 629     if (testIndex_ >= 0 && testIndex_ < static_cast<int>(testPositions_.size()))
 630     {
 631         if (search_.bGrid_)
 632         {
 633             search_.mapPointToGridCell(testPositions_[testIndex], testcell_, xtest_);
 634         }
 635         else
 636         {
 637             copy_rvec(testPositions_[testIndex_], xtest_);
 638         }
 639         if (search_.excls_ != NULL)
 640         {
 641             const int exclIndex  = testExclusionIds_[testIndex];
 642             if (exclIndex < search_.excls_->nr)
 643             {
 644                 const int startIndex = search_.excls_->index[exclIndex];
 645                 nexcl_ = search_.excls_->index[exclIndex + 1] - startIndex;
 646                 excl_  = &search_.excls_->a[startIndex];
 647             }
 648             else
 649             {
 650                 nexcl_ = 0;
 651                 excl_  = NULL;
 652             }
 653         }
 654     }
 655     previ_     = -1;
 656     prevr2_    = 0.0;
 657     exclind_   = 0;
 658     prevnbi_   = 0;
 659     prevcai_   = -1;
 660 }
 661
 662 void AnalysisNeighborhoodPairSearchImpl::nextTestPosition()
 663 {
 664     if (testIndex_ < static_cast<int>(testPositions_.size()))
 665     {
 666         ++testIndex_;
 667         reset(testIndex_);
 668     }
 669 }
 670
 671 bool AnalysisNeighborhoodPairSearchImpl::isExcluded(int j)
 672 {
 673     if (exclind_ < nexcl_)
 674     {
 675         const int refId = search_.refExclusionIds_[j];
 676         while (exclind_ < nexcl_ && excl_[exclind_] < refId)
 677         {
 678             ++exclind_;
 679         }
 680         if (exclind_ < nexcl_ && refId == excl_[exclind_])
 681         {
 682             ++exclind_;
 683             return true;
 684         }
 685     }
 686     return false;
 687 }
 688
 689 void AnalysisNeighborhoodPairSearchImpl::startSearch(
 690         const AnalysisNeighborhoodPositions &positions)
 691 {
 692     testExclusionIds_ = positions.exclusionIds_;
 693     GMX_RELEASE_ASSERT(search_.excls_ == NULL || testExclusionIds_ != NULL,
 694                        "Exclusion IDs must be set when exclusions are enabled");
 695     if (positions.index_ < 0)
 696     {
 697         testPositions_ = constArrayRefFromArray<rvec>(positions.x_, positions.count_);
 698         reset(0);
 699     }
 700     else
 701     {
 702         // Somewhat of a hack: setup the array such that only the last position
 703         // will be used.
 704         testPositions_ = constArrayRefFromArray<rvec>(positions.x_, positions.index_ + 1);
 705         reset(positions.index_);
 706     }
 707 }
 708
 709 template <class Action>
 710 bool AnalysisNeighborhoodPairSearchImpl::searchNext(Action action)
 711 {
 712     while (testIndex_ < static_cast<int>(testPositions_.size()))
 713     {
 714         if (search_.bGrid_)
 715         {
 716             GMX_RELEASE_ASSERT(!search_.bXY_, "Grid-based XY searches not implemented");
 717
 718             int nbi = prevnbi_;
 719             int cai = prevcai_ + 1;
 720
 721             for (; nbi < search_.ngridnb_; ++nbi)
 722             {
 723                 rvec      shift;
 724                 const int ci       = search_.getNeighboringCell(testcell_, nbi, shift);
 725                 const int cellSize = static_cast<int>(search_.cells_[ci].size());
 726                 for (; cai < cellSize; ++cai)
 727                 {
 728                     const int i = search_.cells_[ci][cai];
 729                     if (isExcluded(i))
 730                     {
 731                         continue;
 732                     }
 733                     rvec       dx;
 734                     rvec_sub(xtest_, search_.xref_[i], dx);
 735                     rvec_add(dx, shift, dx);
 736                     const real r2 = norm2(dx);
 737                     if (r2 <= search_.cutoff2_)
 738                     {
 739                         if (action(i, r2))
 740                         {
 741                             prevnbi_ = nbi;
 742                             prevcai_ = cai;
 743                             previ_   = i;
 744                             prevr2_  = r2;
 745                             return true;
 746                         }
 747                     }
 748                 }
 749                 exclind_ = 0;
 750                 cai      = 0;
 751             }
 752         }
 753         else
 754         {
 755             for (int i = previ_ + 1; i < search_.nref_; ++i)
 756             {
 757                 if (isExcluded(i))
 758                 {
 759                     continue;
 760                 }
 761                 rvec dx;
 762                 if (search_.pbc_.ePBC != epbcNONE)
 763                 {
 764                     pbc_dx(&search_.pbc_, xtest_, search_.xref_[i], dx);
 765                 }
 766                 else
 767                 {
 768                     rvec_sub(xtest_, search_.xref_[i], dx);
 769                 }
 770                 const real r2
 771                     = search_.bXY_
 772                         ? dx[XX]*dx[XX] + dx[YY]*dx[YY]
 773                         : norm2(dx);
 774                 if (r2 <= search_.cutoff2_)
 775                 {
 776                     if (action(i, r2))
 777                     {
 778                         previ_  = i;
 779                         prevr2_ = r2;
 780                         return true;
 781                     }
 782                 }
 783             }
 784         }
 785         nextTestPosition();
 786     }
 787     return false;
 788 }
 789
 790 void AnalysisNeighborhoodPairSearchImpl::initFoundPair(
 791         AnalysisNeighborhoodPair *pair) const
 792 {
 793     if (previ_ < 0)
 794     {
 795         *pair = AnalysisNeighborhoodPair();
 796     }
 797     else
 798     {
 799         *pair = AnalysisNeighborhoodPair(previ_, testIndex_, prevr2_);
 800     }
 801 }
 802
 803 }   // namespace internal
 804
 805 namespace
 806 {
 807
 808 /*! \brief
 809  * Search action to find the next neighbor.
 810  *
 811  * Used as the action for AnalysisNeighborhoodPairSearchImpl::searchNext() to
 812  * find the next neighbor.
 813  *
 814  * Simply breaks the loop on the first found neighbor.
 815  */
 816 bool withinAction(int /*i*/, real /*r2*/)
 817 {
 818     return true;
 819 }
 820
 821 /*! \brief
 822  * Search action find the minimum distance.
 823  *
 824  * Used as the action for AnalysisNeighborhoodPairSearchImpl::searchNext() to
 825  * find the nearest neighbor.
 826  *
 827  * With this action, AnalysisNeighborhoodPairSearchImpl::searchNext() always
 828  * returns false, and the output is put into the variables passed by pointer
 829  * into the constructor.  If no neighbors are found, the output variables are
 830  * not modified, i.e., the caller must initialize them.
 831  */
 832 class MindistAction
 833 {
 834     public:
 835         /*! \brief
 836          * Initializes the action with given output locations.
 837          *
 838          * \param[out] closestPoint Index of the closest reference location.
 839          * \param[out] minDist2     Minimum distance squared.
 840          *
 841          * The constructor call does not modify the pointed values, but only
 842          * stores the pointers for later use.
 843          * See the class description for additional semantics.
 844          */
 845         MindistAction(int *closestPoint, real *minDist2)
 846             : closestPoint_(*closestPoint), minDist2_(*minDist2)
 847         {
 848         }
 849
 850         //! Processes a neighbor to find the nearest point.
 851         bool operator()(int i, real r2)
 852         {
 853             if (r2 < minDist2_)
 854             {
 855                 closestPoint_ = i;
 856                 minDist2_     = r2;
 857             }
 858             return false;
 859         }
 860
 861     private:
 862         int     &closestPoint_;
 863         real    &minDist2_;
 864
 865         GMX_DISALLOW_ASSIGN(MindistAction);
 866 };
 867
 868 }   // namespace
 869
 870 /********************************************************************
 871  * AnalysisNeighborhood::Impl
 872  */
 873
 874 class AnalysisNeighborhood::Impl
 875 {
 876     public:
 877         typedef AnalysisNeighborhoodSearch::ImplPointer SearchImplPointer;
 878         typedef std::vector<SearchImplPointer> SearchList;
 879
 880         Impl() : cutoff_(0), excls_(NULL), mode_(eSearchMode_Automatic), bXY_(false)
 881         {
 882         }
 883         ~Impl()
 884         {
 885             SearchList::const_iterator i;
 886             for (i = searchList_.begin(); i != searchList_.end(); ++i)
 887             {
 888                 GMX_RELEASE_ASSERT(i->unique(),
 889                                    "Dangling AnalysisNeighborhoodSearch reference");
 890             }
 891         }
 892
 893         SearchImplPointer getSearch();
 894
 895         tMPI::mutex             createSearchMutex_;
 896         SearchList              searchList_;
 897         real                    cutoff_;
 898         const t_blocka         *excls_;
 899         SearchMode              mode_;
 900         bool                    bXY_;
 901 };
 902
 903 AnalysisNeighborhood::Impl::SearchImplPointer
 904 AnalysisNeighborhood::Impl::getSearch()
 905 {
 906     tMPI::lock_guard<tMPI::mutex> lock(createSearchMutex_);
 907     // TODO: Consider whether this needs to/can be faster, e.g., by keeping a
 908     // separate pool of unused search objects.
 909     SearchList::const_iterator i;
 910     for (i = searchList_.begin(); i != searchList_.end(); ++i)
 911     {
 912         if (i->unique())
 913         {
 914             return *i;
 915         }
 916     }
 917     SearchImplPointer search(new internal::AnalysisNeighborhoodSearchImpl(cutoff_));
 918     searchList_.push_back(search);
 919     return search;
 920 }
 921
 922 /********************************************************************
 923  * AnalysisNeighborhood
 924  */
 925
 926 AnalysisNeighborhood::AnalysisNeighborhood()
 927     : impl_(new Impl)
 928 {
 929 }
 930
 931 AnalysisNeighborhood::~AnalysisNeighborhood()
 932 {
 933 }
 934
 935 void AnalysisNeighborhood::setCutoff(real cutoff)
 936 {
 937     GMX_RELEASE_ASSERT(impl_->searchList_.empty(),
 938                        "Changing the cutoff after initSearch() not currently supported");
 939     impl_->cutoff_ = cutoff;
 940 }
 941
 942 void AnalysisNeighborhood::setXYMode(bool bXY)
 943 {
 944     impl_->bXY_ = bXY;
 945 }
 946
 947 void AnalysisNeighborhood::setTopologyExclusions(const t_blocka *excls)
 948 {
 949     GMX_RELEASE_ASSERT(impl_->searchList_.empty(),
 950                        "Changing the exclusions after initSearch() not currently supported");
 951     impl_->excls_ = excls;
 952 }
 953
 954 void AnalysisNeighborhood::setMode(SearchMode mode)
 955 {
 956     impl_->mode_ = mode;
 957 }
 958
 959 AnalysisNeighborhood::SearchMode AnalysisNeighborhood::mode() const
 960 {
 961     return impl_->mode_;
 962 }
 963
 964 AnalysisNeighborhoodSearch
 965 AnalysisNeighborhood::initSearch(const t_pbc                         *pbc,
 966                                  const AnalysisNeighborhoodPositions &positions)
 967 {
 968     Impl::SearchImplPointer search(impl_->getSearch());
 969     search->init(mode(), impl_->bXY_, impl_->excls_, pbc, positions);
 970     return AnalysisNeighborhoodSearch(search);
 971 }
 972
 973 /********************************************************************
 974  * AnalysisNeighborhoodSearch
 975  */
 976
 977 AnalysisNeighborhoodSearch::AnalysisNeighborhoodSearch()
 978 {
 979 }
 980
 981 AnalysisNeighborhoodSearch::AnalysisNeighborhoodSearch(const ImplPointer &impl)
 982     : impl_(impl)
 983 {
 984 }
 985
 986 void AnalysisNeighborhoodSearch::reset()
 987 {
 988     impl_.reset();
 989 }
 990
 991 AnalysisNeighborhood::SearchMode AnalysisNeighborhoodSearch::mode() const
 992 {
 993     GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
 994     return (impl_->usesGridSearch()
 995             ? AnalysisNeighborhood::eSearchMode_Grid
 996             : AnalysisNeighborhood::eSearchMode_Simple);
 997 }
 998
 999 bool AnalysisNeighborhoodSearch::isWithin(
1000         const AnalysisNeighborhoodPositions &positions) const
1001 {
1002     GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
1003     internal::AnalysisNeighborhoodPairSearchImpl pairSearch(*impl_);
1004     pairSearch.startSearch(positions);
1005     return pairSearch.searchNext(&withinAction);
1006 }
1007
1008 real AnalysisNeighborhoodSearch::minimumDistance(
1009         const AnalysisNeighborhoodPositions &positions) const
1010 {
1011     GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
1012     internal::AnalysisNeighborhoodPairSearchImpl pairSearch(*impl_);
1013     pairSearch.startSearch(positions);
1014     real          minDist2     = impl_->cutoffSquared();
1015     int           closestPoint = -1;
1016     MindistAction action(&closestPoint, &minDist2);
1017     (void)pairSearch.searchNext(action);
1018     return sqrt(minDist2);
1019 }
1020
1021 AnalysisNeighborhoodPair
1022 AnalysisNeighborhoodSearch::nearestPoint(
1023         const AnalysisNeighborhoodPositions &positions) const
1024 {
1025     GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
1026     internal::AnalysisNeighborhoodPairSearchImpl pairSearch(*impl_);
1027     pairSearch.startSearch(positions);
1028     real          minDist2     = impl_->cutoffSquared();
1029     int           closestPoint = -1;
1030     MindistAction action(&closestPoint, &minDist2);
1031     (void)pairSearch.searchNext(action);
1032     return AnalysisNeighborhoodPair(closestPoint, 0, minDist2);
1033 }
1034
1035 AnalysisNeighborhoodPairSearch
1036 AnalysisNeighborhoodSearch::startPairSearch(
1037         const AnalysisNeighborhoodPositions &positions) const
1038 {
1039     GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
1040     Impl::PairSearchImplPointer pairSearch(impl_->getPairSearch());
1041     pairSearch->startSearch(positions);
1042     return AnalysisNeighborhoodPairSearch(pairSearch);
1043 }
1044
1045 /********************************************************************
1046  * AnalysisNeighborhoodPairSearch
1047  */
1048
1049 AnalysisNeighborhoodPairSearch::AnalysisNeighborhoodPairSearch(
1050         const ImplPointer &impl)
1051     : impl_(impl)
1052 {
1053 }
1054
1055 bool AnalysisNeighborhoodPairSearch::findNextPair(AnalysisNeighborhoodPair *pair)
1056 {
1057     bool bFound = impl_->searchNext(&withinAction);
1058     impl_->initFoundPair(pair);
1059     return bFound;
1060 }
1061
1062 void AnalysisNeighborhoodPairSearch::skipRemainingPairsForTestPosition()
1063 {
1064     impl_->nextTestPosition();
1065 }
1066
1067 } // namespace gmx