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