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