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