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