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