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