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