src/gromacs/selection/nbsearch.cpp

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