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