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
  90         explicit AnalysisNeighborhoodSearchImpl(real cutoff);
  91         ~AnalysisNeighborhoodSearchImpl();
  92
  93         void init(AnalysisNeighborhood::SearchMode mode,
  94                   const t_pbc *pbc, int n, const rvec x[]);
  95
  96         //! Calculates offsets to neighboring grid cells that should be considered.
  97         void initGridCellNeighborList();
  98         /*! \brief
  99          * Determines a suitable grid size and sets up the cells.
 100          *
 101          * \param[in]     pbc  Information about the box.
 102          * \returns  false if grid search is not suitable.
 103          */
 104         bool initGridCells(const t_pbc *pbc);
 105         /*! \brief
 106          * Sets ua a search grid for a given box.
 107          *
 108          * \param[in]     pbc  Information about the box.
 109          * \returns  false if grid search is not suitable.
 110          */
 111         bool initGrid(const t_pbc *pbc);
 112         //! Clears all grid cells.
 113         void clearGridCells();
 114         /*! \brief
 115          * Maps a point into a grid cell.
 116          *
 117          * \param[in]  x    Point to map.
 118          * \param[out] cell Indices of the grid cell in which \p x lies.
 119          *
 120          * \p x should be within the triclinic unit cell.
 121          */
 122         void mapPointToGridCell(const rvec x, ivec cell) const;
 123         /*! \brief
 124          * Calculates linear index of a grid cell.
 125          *
 126          * \param[in]  cell Cell indices.
 127          * \returns    Linear index of \p cell.
 128          */
 129         int getGridCellIndex(const ivec cell) const;
 130         /*! \brief
 131          * Adds an index into a grid cell.
 132          *
 133          * \param[in]     cell Cell into which \p i should be added.
 134          * \param[in]     i    Index to add.
 135          */
 136         void addToGridCell(const ivec cell, int i);
 137
 138         /** Whether to try grid searching. */
 139         bool                    bTryGrid_;
 140         /** The cutoff. */
 141         real                    cutoff_;
 142         /** The cutoff squared. */
 143         real                    cutoff2_;
 144
 145         /** Number of reference points for the current frame. */
 146         int                     nref_;
 147         /** Reference point positions. */
 148         const rvec             *xref_;
 149         /** Reference position ids (NULL if not available). */
 150         int                    *refid_;
 151         /** PBC data. */
 152         t_pbc                  *pbc_;
 153
 154         /** Number of excluded reference positions for current test particle. */
 155         int                     nexcl_;
 156         /** Exclusions for current test particle. */
 157         int                    *excl_;
 158
 159         /** Whether grid searching is actually used for the current positions. */
 160         bool                    bGrid_;
 161         /** Array allocated for storing in-unit-cell reference positions. */
 162         rvec                   *xref_alloc_;
 163         /** Allocation count for xref_alloc. */
 164         int                     xref_nalloc_;
 165         /** false if the box is rectangular. */
 166         bool                    bTric_;
 167         /** Box vectors of a single grid cell. */
 168         matrix                  cellbox_;
 169         /** The reciprocal cell vectors as columns; the inverse of \p cellbox. */
 170         matrix                  recipcell_;
 171         /** Number of cells along each dimension. */
 172         ivec                    ncelldim_;
 173         /** Total number of cells. */
 174         int                     ncells_;
 175         /** Number of reference positions in each cell. */
 176         int                    *ncatoms_;
 177         /** List of reference positions in each cell. */
 178         atom_id               **catom_;
 179         /** Allocation counts for each \p catom[i]. */
 180         int                    *catom_nalloc_;
 181         /** Allocation count for the per-cell arrays. */
 182         int                     cells_nalloc_;
 183         /** Number of neighboring cells to consider. */
 184         int                     ngridnb_;
 185         /** Offsets of the neighboring cells to consider. */
 186         ivec                   *gnboffs_;
 187         /** Allocation count for \p gnboffs. */
 188         int                     gnboffs_nalloc_;
 189
 190         PairSearchImplPointer   pairSearch_;
 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         //! Checks whether a reference positiong should be excluded.
 211         bool isExcluded(int j);
 212         //! Searches for the next neighbor.
 213         template <class Action>
 214         bool searchNext(Action action);
 215
 216         const AnalysisNeighborhoodSearchImpl   &search_;
 217         int                                     testIndex_;
 218         /** Stores test position during a pair loop. */
 219         rvec                                    xtest_;
 220         /** Stores the previous returned position during a pair loop. */
 221         int                                     previ_;
 222         /** Stores the current exclusion index during loops. */
 223         int                                     exclind_;
 224         /** Stores the test particle cell index during loops. */
 225         ivec                                    testcell_;
 226         /** Stores the current cell neighbor index during pair loops. */
 227         int                                     prevnbi_;
 228         /** Stores the index within the current cell during pair loops. */
 229         int                                     prevcai_;
 230
 231         GMX_DISALLOW_COPY_AND_ASSIGN(AnalysisNeighborhoodPairSearchImpl);
 232 };
 233
 234 /********************************************************************
 235  * AnalysisNeighborhoodSearchImpl
 236  */
 237
 238 AnalysisNeighborhoodSearchImpl::AnalysisNeighborhoodSearchImpl(real cutoff)
 239     : pairSearch_(new AnalysisNeighborhoodPairSearchImpl(*this))
 240 {
 241     bTryGrid_       = true;
 242     cutoff_         = cutoff;
 243     if (cutoff_ <= 0)
 244     {
 245         cutoff_     = GMX_REAL_MAX;
 246         bTryGrid_   = false;
 247     }
 248     cutoff2_        = sqr(cutoff_);
 249
 250     nref_           = 0;
 251     xref_           = NULL;
 252     refid_          = NULL;
 253     pbc_            = NULL;
 254
 255     nexcl_          = 0;
 256     excl_           = NULL;
 257
 258     bGrid_          = false;
 259
 260     xref_alloc_     = NULL;
 261     xref_nalloc_    = 0;
 262     bTric_          = false;
 263     clear_mat(cellbox_);
 264     clear_mat(recipcell_);
 265     clear_ivec(ncelldim_);
 266     ncells_         = 0;
 267     ncatoms_        = NULL;
 268     catom_          = NULL;
 269     catom_nalloc_   = 0;
 270     cells_nalloc_   = 0;
 271
 272     ngridnb_        = 0;
 273     gnboffs_        = NULL;
 274     gnboffs_nalloc_ = 0;
 275 }
 276
 277 AnalysisNeighborhoodSearchImpl::~AnalysisNeighborhoodSearchImpl()
 278 {
 279     GMX_RELEASE_ASSERT(pairSearch_.unique(),
 280                        "Dangling AnalysisNeighborhoodPairSearch reference");
 281     sfree(xref_alloc_);
 282     sfree(ncatoms_);
 283     if (catom_ != NULL)
 284     {
 285         for (int ci = 0; ci < ncells_; ++ci)
 286         {
 287             sfree(catom_[ci]);
 288         }
 289         sfree(catom_);
 290     }
 291     sfree(catom_nalloc_);
 292     sfree(gnboffs_);
 293 }
 294
 295 void AnalysisNeighborhoodSearchImpl::initGridCellNeighborList()
 296 {
 297     int   maxx, maxy, maxz;
 298     real  rvnorm;
 299
 300     /* Find the extent of the sphere in triclinic coordinates */
 301     maxz   = (int)(cutoff_ * recipcell_[ZZ][ZZ]) + 1;
 302     rvnorm = sqrt(sqr(recipcell_[YY][YY]) + sqr(recipcell_[ZZ][YY]));
 303     maxy   = (int)(cutoff_ * rvnorm) + 1;
 304     rvnorm = sqrt(sqr(recipcell_[XX][XX]) + sqr(recipcell_[YY][XX])
 305                   + sqr(recipcell_[ZZ][XX]));
 306     maxx   = (int)(cutoff_ * rvnorm) + 1;
 307
 308     /* Calculate the number of cells and reallocate if necessary */
 309     ngridnb_ = (2 * maxx + 1) * (2 * maxy + 1) * (2 * maxz + 1);
 310     if (gnboffs_nalloc_ < ngridnb_)
 311     {
 312         gnboffs_nalloc_ = ngridnb_;
 313         srenew(gnboffs_, gnboffs_nalloc_);
 314     }
 315
 316     /* Store the whole cube */
 317     /* TODO: Prune off corners that are not needed */
 318     int i = 0;
 319     for (int x = -maxx; x <= maxx; ++x)
 320     {
 321         for (int y = -maxy; y <= maxy; ++y)
 322         {
 323             for (int z = -maxz; z <= maxz; ++z)
 324             {
 325                 gnboffs_[i][XX] = x;
 326                 gnboffs_[i][YY] = y;
 327                 gnboffs_[i][ZZ] = z;
 328                 ++i;
 329             }
 330         }
 331     }
 332 }
 333
 334 bool AnalysisNeighborhoodSearchImpl::initGridCells(const t_pbc *pbc)
 335 {
 336     const real targetsize =
 337         pow(pbc->box[XX][XX] * pbc->box[YY][YY] * pbc->box[ZZ][ZZ]
 338             * 10 / nref_, static_cast<real>(1./3.));
 339
 340     ncells_ = 1;
 341     for (int dd = 0; dd < DIM; ++dd)
 342     {
 343         ncelldim_[dd] = static_cast<int>(pbc->box[dd][dd] / targetsize);
 344         ncells_      *= ncelldim_[dd];
 345         if (ncelldim_[dd] < 3)
 346         {
 347             return false;
 348         }
 349     }
 350     /* Reallocate if necessary */
 351     if (cells_nalloc_ < ncells_)
 352     {
 353         srenew(ncatoms_, ncells_);
 354         srenew(catom_, ncells_);
 355         srenew(catom_nalloc_, ncells_);
 356         for (int i = cells_nalloc_; i < ncells_; ++i)
 357         {
 358             catom_[i]        = NULL;
 359             catom_nalloc_[i] = 0;
 360         }
 361         cells_nalloc_ = ncells_;
 362     }
 363     return true;
 364 }
 365
 366 bool AnalysisNeighborhoodSearchImpl::initGrid(const t_pbc *pbc)
 367 {
 368     /* TODO: This check could be improved. */
 369     if (0.5*pbc->max_cutoff2 < cutoff2_)
 370     {
 371         return false;
 372     }
 373
 374     if (!initGridCells(pbc))
 375     {
 376         return false;
 377     }
 378
 379     bTric_ = TRICLINIC(pbc->box);
 380     if (bTric_)
 381     {
 382         for (int dd = 0; dd < DIM; ++dd)
 383         {
 384             svmul(1.0 / ncelldim_[dd], pbc->box[dd], cellbox_[dd]);
 385         }
 386         m_inv_ur0(cellbox_, recipcell_);
 387     }
 388     else
 389     {
 390         for (int dd = 0; dd < DIM; ++dd)
 391         {
 392             cellbox_[dd][dd]   = pbc->box[dd][dd] / ncelldim_[dd];
 393             recipcell_[dd][dd] = 1.0 / cellbox_[dd][dd];
 394         }
 395     }
 396     initGridCellNeighborList();
 397     return true;
 398 }
 399
 400 void AnalysisNeighborhoodSearchImpl::mapPointToGridCell(const rvec x,
 401                                                         ivec       cell) const
 402 {
 403     if (bTric_)
 404     {
 405         rvec tx;
 406
 407         tmvmul_ur0(recipcell_, x, tx);
 408         for (int dd = 0; dd < DIM; ++dd)
 409         {
 410             cell[dd] = static_cast<int>(tx[dd]);
 411         }
 412     }
 413     else
 414     {
 415         for (int dd = 0; dd < DIM; ++dd)
 416         {
 417             cell[dd] = static_cast<int>(x[dd] * recipcell_[dd][dd]);
 418         }
 419     }
 420 }
 421
 422 int AnalysisNeighborhoodSearchImpl::getGridCellIndex(const ivec cell) const
 423 {
 424     return cell[XX] + cell[YY] * ncelldim_[XX]
 425            + cell[ZZ] * ncelldim_[XX] * ncelldim_[YY];
 426 }
 427
 428 void AnalysisNeighborhoodSearchImpl::clearGridCells()
 429 {
 430     for (int ci = 0; ci < ncells_; ++ci)
 431     {
 432         ncatoms_[ci] = 0;
 433     }
 434 }
 435
 436 void AnalysisNeighborhoodSearchImpl::addToGridCell(const ivec cell, int i)
 437 {
 438     const int ci = getGridCellIndex(cell);
 439
 440     if (ncatoms_[ci] == catom_nalloc_[ci])
 441     {
 442         catom_nalloc_[ci] += 10;
 443         srenew(catom_[ci], catom_nalloc_[ci]);
 444     }
 445     catom_[ci][ncatoms_[ci]++] = i;
 446 }
 447
 448 void AnalysisNeighborhoodSearchImpl::init(
 449         AnalysisNeighborhood::SearchMode mode,
 450         const t_pbc *pbc, int n, const rvec x[])
 451 {
 452     pbc_  = const_cast<t_pbc *>(pbc);
 453     nref_ = n;
 454     // TODO: Consider whether it would be possible to support grid searching in
 455     // more cases.
 456     if (mode == AnalysisNeighborhood::eSearchMode_Simple
 457         || pbc_ == NULL || pbc_->ePBC != epbcXYZ)
 458     {
 459         bGrid_ = false;
 460     }
 461     else if (bTryGrid_)
 462     {
 463         // TODO: Actually implement forcing eSearchMode_Grid
 464         bGrid_ = initGrid(pbc_);
 465     }
 466     if (bGrid_)
 467     {
 468         if (xref_nalloc_ < n)
 469         {
 470             srenew(xref_alloc_, n);
 471             xref_nalloc_ = n;
 472         }
 473         xref_ = xref_alloc_;
 474         clearGridCells();
 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_ = NULL;
 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             /* TODO: Calculate the required PBC shift outside the inner loop */
 595             for (; cai < search_.ncatoms_[ci]; ++cai)
 596             {
 597                 const int i = search_.catom_[ci][cai];
 598                 if (isExcluded(i))
 599                 {
 600                     continue;
 601                 }
 602                 rvec       dx;
 603                 pbc_dx_aiuc(search_.pbc_, xtest_, search_.xref_[i], dx);
 604                 const real r2 = norm2(dx);
 605                 if (r2 <= search_.cutoff2_)
 606                 {
 607                     if (action(i, r2))
 608                     {
 609                         prevnbi_ = nbi;
 610                         prevcai_ = cai;
 611                         previ_   = i;
 612                         return true;
 613                     }
 614                 }
 615             }
 616             exclind_ = 0;
 617             cai      = 0;
 618         }
 619     }
 620     else
 621     {
 622         for (int i = previ_ + 1; i < search_.nref_; ++i)
 623         {
 624             if (isExcluded(i))
 625             {
 626                 continue;
 627             }
 628             rvec dx;
 629             if (search_.pbc_)
 630             {
 631                 pbc_dx(search_.pbc_, xtest_, search_.xref_[i], dx);
 632             }
 633             else
 634             {
 635                 rvec_sub(xtest_, search_.xref_[i], dx);
 636             }
 637             const real r2 = norm2(dx);
 638             if (r2 <= search_.cutoff2_)
 639             {
 640                 if (action(i, r2))
 641                 {
 642                     previ_ = i;
 643                     return true;
 644                 }
 645             }
 646         }
 647     }
 648     return false;
 649 }
 650
 651 }   // namespace internal
 652
 653 namespace
 654 {
 655
 656 /*! \brief
 657  * Search action to find the next neighbor.
 658  *
 659  * Used as the action for AnalysisNeighborhoodPairSearchImpl::searchNext() to
 660  * find the next neighbor.
 661  *
 662  * Simply breaks the loop on the first found neighbor.
 663  */
 664 bool withinAction(int /*i*/, real /*r2*/)
 665 {
 666     return true;
 667 }
 668
 669 /*! \brief
 670  * Search action find the minimum distance.
 671  *
 672  * Used as the action for AnalysisNeighborhoodPairSearchImpl::searchNext() to
 673  * find the nearest neighbor.
 674  *
 675  * With this action, AnalysisNeighborhoodPairSearchImpl::searchNext() always
 676  * returns false, and the output is put into the variables passed by pointer
 677  * into the constructor.  If no neighbors are found, the output variables are
 678  * not modified, i.e., the caller must initialize them.
 679  */
 680 class MindistAction
 681 {
 682     public:
 683         /*! \brief
 684          * Initializes the action with given output locations.
 685          *
 686          * \param[out] closestPoint Index of the closest reference location.
 687          * \param[out] minDist2     Minimum distance squared.
 688          *
 689          * The constructor call does not modify the pointed values, but only
 690          * stores the pointers for later use.
 691          * See the class description for additional semantics.
 692          */
 693         MindistAction(int *closestPoint, real *minDist2)
 694             : closestPoint_(*closestPoint), minDist2_(*minDist2)
 695         {
 696         }
 697
 698         //! Processes a neighbor to find the nearest point.
 699         bool operator()(int i, real r2)
 700         {
 701             if (r2 < minDist2_)
 702             {
 703                 closestPoint_ = i;
 704                 minDist2_     = r2;
 705             }
 706             return false;
 707         }
 708
 709     private:
 710         int     &closestPoint_;
 711         real    &minDist2_;
 712
 713         GMX_DISALLOW_ASSIGN(MindistAction);
 714 };
 715
 716 }   // namespace
 717
 718 /********************************************************************
 719  * AnalysisNeighborhood::Impl
 720  */
 721
 722 class AnalysisNeighborhood::Impl
 723 {
 724     public:
 725         typedef AnalysisNeighborhoodSearch::ImplPointer SearchImplPointer;
 726         typedef std::vector<SearchImplPointer> SearchList;
 727
 728         Impl() : cutoff_(0), mode_(eSearchMode_Automatic)
 729         {
 730         }
 731         ~Impl()
 732         {
 733             SearchList::const_iterator i;
 734             for (i = searchList_.begin(); i != searchList_.end(); ++i)
 735             {
 736                 GMX_RELEASE_ASSERT(i->unique(),
 737                                    "Dangling AnalysisNeighborhoodSearch reference");
 738             }
 739         }
 740
 741         SearchImplPointer getSearch();
 742
 743         tMPI::mutex             createSearchMutex_;
 744         SearchList              searchList_;
 745         real                    cutoff_;
 746         SearchMode              mode_;
 747 };
 748
 749 AnalysisNeighborhood::Impl::SearchImplPointer
 750 AnalysisNeighborhood::Impl::getSearch()
 751 {
 752     tMPI::lock_guard<tMPI::mutex> lock(createSearchMutex_);
 753     // TODO: Consider whether this needs to/can be faster, e.g., by keeping a
 754     // separate pool of unused search objects.
 755     SearchList::const_iterator i;
 756     for (i = searchList_.begin(); i != searchList_.end(); ++i)
 757     {
 758         if (i->unique())
 759         {
 760             return *i;
 761         }
 762     }
 763     SearchImplPointer search(new internal::AnalysisNeighborhoodSearchImpl(cutoff_));
 764     searchList_.push_back(search);
 765     return search;
 766 }
 767
 768 /********************************************************************
 769  * AnalysisNeighborhood
 770  */
 771
 772 AnalysisNeighborhood::AnalysisNeighborhood()
 773     : impl_(new Impl)
 774 {
 775 }
 776
 777 AnalysisNeighborhood::~AnalysisNeighborhood()
 778 {
 779 }
 780
 781 void AnalysisNeighborhood::setCutoff(real cutoff)
 782 {
 783     GMX_RELEASE_ASSERT(impl_->searchList_.empty(),
 784                        "Changing the cutoff after initSearch() not currently supported");
 785     impl_->cutoff_ = cutoff;
 786 }
 787
 788 void AnalysisNeighborhood::setMode(SearchMode mode)
 789 {
 790     impl_->mode_ = mode;
 791 }
 792
 793 AnalysisNeighborhood::SearchMode AnalysisNeighborhood::mode() const
 794 {
 795     return impl_->mode_;
 796 }
 797
 798 AnalysisNeighborhoodSearch
 799 AnalysisNeighborhood::initSearch(const t_pbc *pbc, int n, const rvec x[])
 800 {
 801     Impl::SearchImplPointer search(impl_->getSearch());
 802     search->init(mode(), pbc, n, x);
 803     return AnalysisNeighborhoodSearch(search);
 804 }
 805
 806 AnalysisNeighborhoodSearch
 807 AnalysisNeighborhood::initSearch(const t_pbc *pbc, const gmx_ana_pos_t *p)
 808 {
 809     Impl::SearchImplPointer search(impl_->getSearch());
 810     search->init(mode(), pbc, p->nr, p->x);
 811     search->refid_ = p->m.refid;
 812     return AnalysisNeighborhoodSearch(search);
 813 }
 814
 815 /********************************************************************
 816  * AnalysisNeighborhoodSearch
 817  */
 818
 819 AnalysisNeighborhoodSearch::AnalysisNeighborhoodSearch()
 820 {
 821 }
 822
 823 AnalysisNeighborhoodSearch::AnalysisNeighborhoodSearch(const ImplPointer &impl)
 824     : impl_(impl)
 825 {
 826 }
 827
 828 void AnalysisNeighborhoodSearch::reset()
 829 {
 830     impl_.reset();
 831 }
 832
 833 AnalysisNeighborhood::SearchMode AnalysisNeighborhoodSearch::mode() const
 834 {
 835     GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
 836     return (impl_->bGrid_
 837             ? AnalysisNeighborhood::eSearchMode_Grid
 838             : AnalysisNeighborhood::eSearchMode_Simple);
 839 }
 840
 841 bool AnalysisNeighborhoodSearch::isWithin(const rvec x) const
 842 {
 843     GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
 844     internal::AnalysisNeighborhoodPairSearchImpl pairSearch(*impl_);
 845     pairSearch.startSearch(x, 0);
 846     return pairSearch.searchNext(&withinAction);
 847 }
 848
 849 bool AnalysisNeighborhoodSearch::isWithin(const gmx_ana_pos_t *p, int i) const
 850 {
 851     return isWithin(p->x[i]);
 852 }
 853
 854 real AnalysisNeighborhoodSearch::minimumDistance(const rvec x) const
 855 {
 856     GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
 857     internal::AnalysisNeighborhoodPairSearchImpl pairSearch(*impl_);
 858     pairSearch.startSearch(x, 0);
 859     real          minDist2     = impl_->cutoff2_;
 860     int           closestPoint = -1;
 861     MindistAction action(&closestPoint, &minDist2);
 862     (void)pairSearch.searchNext(action);
 863     return sqrt(minDist2);
 864 }
 865
 866 real AnalysisNeighborhoodSearch::minimumDistance(const gmx_ana_pos_t *p, int i) const
 867 {
 868     return minimumDistance(p->x[i]);
 869 }
 870
 871 AnalysisNeighborhoodPair
 872 AnalysisNeighborhoodSearch::nearestPoint(const rvec x) const
 873 {
 874     GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
 875     internal::AnalysisNeighborhoodPairSearchImpl pairSearch(*impl_);
 876     pairSearch.startSearch(x, 0);
 877     real          minDist2     = impl_->cutoff2_;
 878     int           closestPoint = -1;
 879     MindistAction action(&closestPoint, &minDist2);
 880     (void)pairSearch.searchNext(action);
 881     return AnalysisNeighborhoodPair(closestPoint, 0);
 882 }
 883
 884 AnalysisNeighborhoodPair
 885 AnalysisNeighborhoodSearch::nearestPoint(const gmx_ana_pos_t *p, int i) const
 886 {
 887     GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
 888     internal::AnalysisNeighborhoodPairSearchImpl pairSearch(*impl_);
 889     pairSearch.startSearch(p->x[i], 0);
 890     real          minDist2     = impl_->cutoff2_;
 891     int           closestPoint = -1;
 892     MindistAction action(&closestPoint, &minDist2);
 893     (void)pairSearch.searchNext(action);
 894     return AnalysisNeighborhoodPair(closestPoint, i);
 895 }
 896
 897 AnalysisNeighborhoodPairSearch
 898 AnalysisNeighborhoodSearch::startPairSearch(const rvec x)
 899 {
 900     GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
 901     GMX_RELEASE_ASSERT(impl_->pairSearch_.unique(),
 902                        "Multiple concurrent searches not currently supported");
 903     impl_->pairSearch_->startSearch(x, 0);
 904     return AnalysisNeighborhoodPairSearch(impl_->pairSearch_);
 905 }
 906
 907 AnalysisNeighborhoodPairSearch
 908 AnalysisNeighborhoodSearch::startPairSearch(const gmx_ana_pos_t *p, int i)
 909 {
 910     GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
 911     GMX_RELEASE_ASSERT(impl_->pairSearch_.unique(),
 912                        "Multiple concurrent searches not currently supported");
 913     impl_->pairSearch_->startSearch(p->x[i], i);
 914     return AnalysisNeighborhoodPairSearch(impl_->pairSearch_);
 915 }
 916
 917 /********************************************************************
 918  * AnalysisNeighborhoodPairSearch
 919  */
 920
 921 AnalysisNeighborhoodPairSearch::AnalysisNeighborhoodPairSearch(
 922         const ImplPointer &impl)
 923     : impl_(impl)
 924 {
 925 }
 926
 927 bool AnalysisNeighborhoodPairSearch::findNextPair(AnalysisNeighborhoodPair *pair)
 928 {
 929     if (impl_->searchNext(&withinAction))
 930     {
 931         *pair = AnalysisNeighborhoodPair(impl_->previ_, impl_->testIndex_);
 932         return true;
 933     }
 934     *pair = AnalysisNeighborhoodPair();
 935     return false;
 936 }
 937
 938 } // namespace gmx