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 /*! \page nbsearch Neighborhood search routines
  36  *
  37  * Functions to find particles within a neighborhood of a set of particles
  38  * are defined in nbsearch.h.
  39  * The usage is simple: a data structure is allocated with
  40  * gmx_ana_nbsearch_create(), and the box shape and reference positions for a
  41  * frame are set using gmx_ana_nbsearch_init() or gmx_ana_nbsearch_pos_init().
  42  * Searches can then be performed with gmx_ana_nbsearch_is_within() and
  43  * gmx_ana_nbsearch_mindist(), or with versions that take the \c gmx_ana_pos_t
  44  * data structure.
  45  * When the data structure is no longer required, it can be freed with
  46  * gmx_ana_nbsearch_free().
  47  *
  48  * \internal
  49  *
  50  * \todo
  51  * The grid implementation could still be optimized in several different ways:
  52  *   - Triclinic grid cells are not the most efficient shape, but make PBC
  53  *     handling easier.
  54  *   - Precalculating the required PBC shift for a pair of cells outside the
  55  *     inner loop. After this is done, it should be quite straightforward to
  56  *     move to rectangular cells.
  57  *   - Pruning grid cells from the search list if they are completely outside
  58  *     the sphere that is being considered.
  59  *   - A better heuristic could be added for falling back to simple loops for a
  60  *     small number of reference particles.
  61  *   - A better heuristic for selecting the grid size.
  62  *   - A multi-level grid implementation could be used to be able to use small
  63  *     grids for short cutoffs with very inhomogeneous particle distributions
  64  *     without a memory cost.
  65  */
  66 /*! \internal \file
  67  * \brief
  68  * Implements functions in nbsearch.h.
  69  *
  70  * \author Teemu Murtola <teemu.murtola@gmail.com>
  71  * \ingroup module_selection
  72  */
  73 #ifdef HAVE_CONFIG_H
  74 #include "config.h"
  75 #endif
  76
  77 #include <math.h>
  78
  79 #include <algorithm>
  80
  81 #include "gromacs/legacyheaders/smalloc.h"
  82 #include "gromacs/legacyheaders/typedefs.h"
  83 #include "gromacs/legacyheaders/pbc.h"
  84 #include "gromacs/legacyheaders/vec.h"
  85
  86 #include "gromacs/selection/nbsearch.h"
  87 #include "gromacs/selection/position.h"
  88 #include "gromacs/utility/gmxassert.h"
  89
  90 /*! \internal \brief
  91  * Data structure for neighborhood searches.
  92  */
  93 struct gmx_ana_nbsearch_t
  94 {
  95     /** The cutoff. */
  96     real           cutoff;
  97     /** The cutoff squared. */
  98     real           cutoff2;
  99     /** Maximum number of reference points. */
 100     int            maxnref;
 101
 102     /** Number of reference points for the current frame. */
 103     int            nref;
 104     /** Reference point positions. */
 105     rvec          *xref;
 106     /** Reference position ids (NULL if not available). */
 107     int           *refid;
 108     /** PBC data. */
 109     t_pbc         *pbc;
 110
 111     /** Number of excluded reference positions for current test particle. */
 112     int            nexcl;
 113     /** Exclusions for current test particle. */
 114     int           *excl;
 115
 116     /** Whether to try grid searching. */
 117     bool           bTryGrid;
 118     /** Whether grid searching is actually used for the current positions. */
 119     bool           bGrid;
 120     /** Array allocated for storing in-unit-cell reference positions. */
 121     rvec          *xref_alloc;
 122     /** Allocation count for xref_alloc. */
 123     int            xref_nalloc;
 124     /** false if the box is rectangular. */
 125     bool           bTric;
 126     /** Box vectors of a single grid cell. */
 127     matrix         cellbox;
 128     /** The reciprocal cell vectors as columns; the inverse of \p cellbox. */
 129     matrix         recipcell;
 130     /** Number of cells along each dimension. */
 131     ivec           ncelldim;
 132     /** Total number of cells. */
 133     int            ncells;
 134     /** Number of reference positions in each cell. */
 135     int           *ncatoms;
 136     /** List of reference positions in each cell. */
 137     atom_id      **catom;
 138     /** Allocation counts for each \p catom[i]. */
 139     int           *catom_nalloc;
 140     /** Allocation count for the per-cell arrays. */
 141     int            cells_nalloc;
 142     /** Number of neighboring cells to consider. */
 143     int            ngridnb;
 144     /** Offsets of the neighboring cells to consider. */
 145     ivec          *gnboffs;
 146     /** Allocation count for \p gnboffs. */
 147     int            gnboffs_nalloc;
 148
 149     /** Stores test position during a pair loop. */
 150     rvec           xtest;
 151     /** Stores the previous returned position during a pair loop. */
 152     int            previ;
 153     /** Stores the current exclusion index during loops. */
 154     int            exclind;
 155     /** Stores the test particle cell index during loops. */
 156     ivec           testcell;
 157     /** Stores the current cell neighbor index during pair loops. */
 158     int            prevnbi;
 159     /** Stores the index within the current cell during pair loops. */
 160     int            prevcai;
 161 };
 162
 163 /*!
 164  * \param[in]  cutoff Cutoff distance for the search
 165  *   (<=0 stands for no cutoff).
 166  * \param[in]  maxn   Maximum number of reference particles.
 167  * \returns  Created neighborhood search data structure.
 168  */
 169 gmx_ana_nbsearch_t *
 170 gmx_ana_nbsearch_create(real cutoff, int maxn)
 171 {
 172     gmx_ana_nbsearch_t *d;
 173
 174     snew(d, 1);
 175     d->bTryGrid = true;
 176     if (cutoff <= 0)
 177     {
 178         cutoff      = GMX_REAL_MAX;
 179         d->bTryGrid = false;
 180     }
 181     d->cutoff  = cutoff;
 182     d->cutoff2 = sqr(cutoff);
 183     d->maxnref = maxn;
 184
 185     d->xref    = NULL;
 186     d->nexcl   = 0;
 187     d->exclind = 0;
 188
 189     d->xref_alloc   = NULL;
 190     d->xref_nalloc  = 0;
 191     d->ncells       = 0;
 192     d->ncatoms      = NULL;
 193     d->catom        = NULL;
 194     d->catom_nalloc = 0;
 195     d->cells_nalloc = 0;
 196
 197     d->ngridnb        = 0;
 198     d->gnboffs        = NULL;
 199     d->gnboffs_nalloc = 0;
 200
 201     return d;
 202 }
 203
 204 /*!
 205  * \param     d Data structure to free.
 206  *
 207  * After the call, the pointer \p d is no longer valid.
 208  */
 209 void
 210 gmx_ana_nbsearch_free(gmx_ana_nbsearch_t *d)
 211 {
 212     sfree(d->xref_alloc);
 213     sfree(d->ncatoms);
 214     if (d->catom)
 215     {
 216         int ci;
 217
 218         for (ci = 0; ci < d->ncells; ++ci)
 219         {
 220             sfree(d->catom[ci]);
 221         }
 222         sfree(d->catom);
 223     }
 224     sfree(d->catom_nalloc);
 225     sfree(d->gnboffs);
 226     sfree(d);
 227 }
 228
 229 /*! \brief
 230  * Calculates offsets to neighboring grid cells that should be considered.
 231  *
 232  * \param[in,out] d    Grid information.
 233  * \param[in]     pbc  Information about the box.
 234  */
 235 static void
 236 grid_init_cell_nblist(gmx_ana_nbsearch_t *d, t_pbc *pbc)
 237 {
 238     int   maxx, maxy, maxz;
 239     int   x, y, z, i;
 240     real  rvnorm;
 241
 242     /* Find the extent of the sphere in triclinic coordinates */
 243     maxz   = (int)(d->cutoff * d->recipcell[ZZ][ZZ]) + 1;
 244     rvnorm = sqrt(sqr(d->recipcell[YY][YY]) + sqr(d->recipcell[ZZ][YY]));
 245     maxy   = (int)(d->cutoff * rvnorm) + 1;
 246     rvnorm = sqrt(sqr(d->recipcell[XX][XX]) + sqr(d->recipcell[YY][XX])
 247                   + sqr(d->recipcell[ZZ][XX]));
 248     maxx = (int)(d->cutoff * rvnorm) + 1;
 249
 250     /* Calculate the number of cells and reallocate if necessary */
 251     d->ngridnb = (2 * maxx + 1) * (2 * maxy + 1) * (2 * maxz + 1);
 252     if (d->gnboffs_nalloc < d->ngridnb)
 253     {
 254         d->gnboffs_nalloc = d->ngridnb;
 255         srenew(d->gnboffs, d->gnboffs_nalloc);
 256     }
 257
 258     /* Store the whole cube */
 259     /* TODO: Prune off corners that are not needed */
 260     i = 0;
 261     for (x = -maxx; x <= maxx; ++x)
 262     {
 263         for (y = -maxy; y <= maxy; ++y)
 264         {
 265             for (z = -maxz; z <= maxz; ++z)
 266             {
 267                 d->gnboffs[i][XX] = x;
 268                 d->gnboffs[i][YY] = y;
 269                 d->gnboffs[i][ZZ] = z;
 270                 ++i;
 271             }
 272         }
 273     }
 274 }
 275
 276 /*! \brief
 277  * Determines a suitable grid size.
 278  *
 279  * \param[in,out] d    Grid information.
 280  * \param[in]     pbc  Information about the box.
 281  * \returns  false if grid search is not suitable.
 282  */
 283 static bool
 284 grid_setup_cells(gmx_ana_nbsearch_t *d, t_pbc *pbc)
 285 {
 286     real targetsize;
 287     int  dd;
 288
 289 #ifdef HAVE_CBRT
 290     targetsize = cbrt(pbc->box[XX][XX] * pbc->box[YY][YY] * pbc->box[ZZ][ZZ]
 291                       * 10 / d->nref);
 292 #else
 293     targetsize = pow(pbc->box[XX][XX] * pbc->box[YY][YY] * pbc->box[ZZ][ZZ]
 294                      * 10 / d->nref, (real)(1./3.));
 295 #endif
 296
 297     d->ncells = 1;
 298     for (dd = 0; dd < DIM; ++dd)
 299     {
 300         d->ncelldim[dd] = (int)(pbc->box[dd][dd] / targetsize);
 301         d->ncells      *= d->ncelldim[dd];
 302         if (d->ncelldim[dd] < 3)
 303         {
 304             return false;
 305         }
 306     }
 307     /* Reallocate if necessary */
 308     if (d->cells_nalloc < d->ncells)
 309     {
 310         int  i;
 311
 312         srenew(d->ncatoms, d->ncells);
 313         srenew(d->catom, d->ncells);
 314         srenew(d->catom_nalloc, d->ncells);
 315         for (i = d->cells_nalloc; i < d->ncells; ++i)
 316         {
 317             d->catom[i]        = NULL;
 318             d->catom_nalloc[i] = 0;
 319         }
 320         d->cells_nalloc = d->ncells;
 321     }
 322     return true;
 323 }
 324
 325 /*! \brief
 326  * Sets ua a search grid for a given box.
 327  *
 328  * \param[in,out] d    Grid information.
 329  * \param[in]     pbc  Information about the box.
 330  * \returns  false if grid search is not suitable.
 331  */
 332 static bool
 333 grid_set_box(gmx_ana_nbsearch_t *d, t_pbc *pbc)
 334 {
 335     int dd;
 336
 337     /* TODO: This check could be improved. */
 338     if (0.5*pbc->max_cutoff2 < d->cutoff2)
 339     {
 340         return false;
 341     }
 342
 343     if (!grid_setup_cells(d, pbc))
 344     {
 345         return false;
 346     }
 347
 348     d->bTric = TRICLINIC(pbc->box);
 349     if (d->bTric)
 350     {
 351         for (dd = 0; dd < DIM; ++dd)
 352         {
 353             svmul(1.0 / d->ncelldim[dd], pbc->box[dd], d->cellbox[dd]);
 354         }
 355         m_inv_ur0(d->cellbox, d->recipcell);
 356     }
 357     else
 358     {
 359         for (dd = 0; dd < DIM; ++dd)
 360         {
 361             d->cellbox[dd][dd]   = pbc->box[dd][dd] / d->ncelldim[dd];
 362             d->recipcell[dd][dd] = 1 / d->cellbox[dd][dd];
 363         }
 364     }
 365     grid_init_cell_nblist(d, pbc);
 366     return true;
 367 }
 368
 369 /*! \brief
 370  * Maps a point into a grid cell.
 371  *
 372  * \param[in]  d    Grid information.
 373  * \param[in]  x    Point to map.
 374  * \param[out] cell Indices of the grid cell in which \p x lies.
 375  *
 376  * \p x should be in the triclinic unit cell.
 377  */
 378 static void
 379 grid_map_onto(gmx_ana_nbsearch_t *d, const rvec x, ivec cell)
 380 {
 381     int dd;
 382
 383     if (d->bTric)
 384     {
 385         rvec tx;
 386
 387         tmvmul_ur0(d->recipcell, x, tx);
 388         for (dd = 0; dd < DIM; ++dd)
 389         {
 390             cell[dd] = (int)tx[dd];
 391         }
 392     }
 393     else
 394     {
 395         for (dd = 0; dd < DIM; ++dd)
 396         {
 397             cell[dd] = (int)(x[dd] * d->recipcell[dd][dd]);
 398         }
 399     }
 400 }
 401
 402 /*! \brief
 403  * Calculates linear index of a grid cell.
 404  *
 405  * \param[in]  d    Grid information.
 406  * \param[in]  cell Cell indices.
 407  * \returns    Linear index of \p cell.
 408  */
 409 static int
 410 grid_index(gmx_ana_nbsearch_t *d, const ivec cell)
 411 {
 412     return cell[XX] + cell[YY] * d->ncelldim[XX]
 413            + cell[ZZ] * d->ncelldim[XX] * d->ncelldim[YY];
 414 }
 415
 416 /*! \brief
 417  * Clears all grid cells.
 418  *
 419  * \param[in,out] d    Grid information.
 420  */
 421 static void
 422 grid_clear_cells(gmx_ana_nbsearch_t *d)
 423 {
 424     int  ci;
 425
 426     for (ci = 0; ci < d->ncells; ++ci)
 427     {
 428         d->ncatoms[ci] = 0;
 429     }
 430 }
 431
 432 /*! \brief
 433  * Adds an index into a grid cell.
 434  *
 435  * \param[in,out] d    Grid information.
 436  * \param[in]     cell Cell into which \p i should be added.
 437  * \param[in]     i    Index to add.
 438  */
 439 static void
 440 grid_add_to_cell(gmx_ana_nbsearch_t *d, const ivec cell, int i)
 441 {
 442     int ci = grid_index(d, cell);
 443
 444     if (d->ncatoms[ci] == d->catom_nalloc[ci])
 445     {
 446         d->catom_nalloc[ci] += 10;
 447         srenew(d->catom[ci], d->catom_nalloc[ci]);
 448     }
 449     d->catom[ci][d->ncatoms[ci]++] = i;
 450 }
 451
 452 /*!
 453  * \param[in,out] d   Neighborhood search data structure.
 454  * \param[in]     pbc PBC information for the frame.
 455  * \param[in]     n   Number of reference positions for the frame.
 456  * \param[in]     x   \p n reference positions for the frame.
 457  *
 458  * Initializes the data structure \p d such that it can be used to search
 459  * for the neighbors of \p x.
 460  */
 461 void
 462 gmx_ana_nbsearch_init(gmx_ana_nbsearch_t *d, t_pbc *pbc, int n, const rvec x[])
 463 {
 464     d->pbc  = pbc;
 465     d->nref = n;
 466     // TODO: Consider whether it would be possible to support grid searching in
 467     // more cases.
 468     if (pbc == NULL || pbc->ePBC != epbcXYZ)
 469     {
 470         d->bGrid = false;
 471     }
 472     else if (d->bTryGrid)
 473     {
 474         d->bGrid = grid_set_box(d, pbc);
 475     }
 476     if (d->bGrid)
 477     {
 478         int  i;
 479
 480         if (d->xref_nalloc < n)
 481         {
 482             const int allocCount = std::max(d->maxnref, n);
 483             srenew(d->xref_alloc, allocCount);
 484             d->xref_nalloc = allocCount;
 485         }
 486         d->xref = d->xref_alloc;
 487         grid_clear_cells(d);
 488
 489         for (i = 0; i < n; ++i)
 490         {
 491             copy_rvec(x[i], d->xref[i]);
 492         }
 493         put_atoms_in_triclinic_unitcell(ecenterTRIC, pbc->box, n, d->xref);
 494         for (i = 0; i < n; ++i)
 495         {
 496             ivec refcell;
 497
 498             grid_map_onto(d, d->xref[i], refcell);
 499             grid_add_to_cell(d, refcell, i);
 500         }
 501     }
 502     else
 503     {
 504         // Won't be modified in this case, but when a grid is used,
 505         // xref _is_ modified, so it can't be const.
 506         d->xref = const_cast<rvec *>(x);
 507     }
 508     d->refid = NULL;
 509 }
 510
 511 /*!
 512  * \param[in,out] d   Neighborhood search data structure.
 513  * \param[in]     pbc PBC information for the frame.
 514  * \param[in]     p   Reference positions for the frame.
 515  *
 516  * A convenience wrapper for gmx_ana_nbsearch_init().
 517  */
 518 void
 519 gmx_ana_nbsearch_pos_init(gmx_ana_nbsearch_t *d, t_pbc *pbc, const gmx_ana_pos_t *p)
 520 {
 521     gmx_ana_nbsearch_init(d, pbc, p->nr, p->x);
 522     d->refid = p->m.refid;
 523 }
 524
 525 /*!
 526  * \param[in,out] d     Neighborhood search data structure.
 527  * \param[in]     nexcl Number of reference positions to exclude from next
 528  *      search.
 529  * \param[in]     excl  Indices of reference positions to exclude.
 530  *
 531  * The set exclusions remain in effect until the next call of this function.
 532  */
 533 void
 534 gmx_ana_nbsearch_set_excl(gmx_ana_nbsearch_t *d, int nexcl, int excl[])
 535 {
 536
 537     d->nexcl = nexcl;
 538     d->excl  = excl;
 539 }
 540
 541 /*! \brief
 542  * Helper function to check whether a reference point should be excluded.
 543  */
 544 static bool
 545 is_excluded(gmx_ana_nbsearch_t *d, int j)
 546 {
 547     if (d->exclind < d->nexcl)
 548     {
 549         if (d->refid)
 550         {
 551             while (d->exclind < d->nexcl &&d->refid[j] > d->excl[d->exclind])
 552             {
 553                 ++d->exclind;
 554             }
 555             if (d->exclind < d->nexcl && d->refid[j] == d->excl[d->exclind])
 556             {
 557                 ++d->exclind;
 558                 return true;
 559             }
 560         }
 561         else
 562         {
 563             while (d->bGrid && d->exclind < d->nexcl && d->excl[d->exclind] < j)
 564             {
 565                 ++d->exclind;
 566             }
 567             if (d->excl[d->exclind] == j)
 568             {
 569                 ++d->exclind;
 570                 return true;
 571             }
 572         }
 573     }
 574     return false;
 575 }
 576
 577 /*! \brief
 578  * Initializes a grid search to find reference positions neighboring \p x.
 579  */
 580 static void
 581 grid_search_start(gmx_ana_nbsearch_t *d, const rvec x)
 582 {
 583     copy_rvec(x, d->xtest);
 584     if (d->bGrid)
 585     {
 586         put_atoms_in_triclinic_unitcell(ecenterTRIC, d->pbc->box, 1, &d->xtest);
 587         grid_map_onto(d, d->xtest, d->testcell);
 588         d->prevnbi = 0;
 589         d->prevcai = -1;
 590     }
 591     d->previ   = -1;
 592     d->exclind = 0;
 593 }
 594
 595 /*! \brief
 596  * Does a grid search.
 597  */
 598 static bool
 599 grid_search(gmx_ana_nbsearch_t *d,
 600             bool (*action)(gmx_ana_nbsearch_t *d, int i, real r2))
 601 {
 602     int  i;
 603     rvec dx;
 604     real r2;
 605
 606     if (d->bGrid)
 607     {
 608         int  nbi, ci, cai;
 609
 610         nbi = d->prevnbi;
 611         cai = d->prevcai + 1;
 612
 613         for (; nbi < d->ngridnb; ++nbi)
 614         {
 615             ivec cell;
 616
 617             ivec_add(d->testcell, d->gnboffs[nbi], cell);
 618             cell[XX] = (cell[XX] + d->ncelldim[XX]) % d->ncelldim[XX];
 619             cell[YY] = (cell[YY] + d->ncelldim[YY]) % d->ncelldim[YY];
 620             cell[ZZ] = (cell[ZZ] + d->ncelldim[ZZ]) % d->ncelldim[ZZ];
 621             ci       = grid_index(d, cell);
 622             /* TODO: Calculate the required PBC shift outside the inner loop */
 623             for (; cai < d->ncatoms[ci]; ++cai)
 624             {
 625                 i = d->catom[ci][cai];
 626                 if (is_excluded(d, i))
 627                 {
 628                     continue;
 629                 }
 630                 pbc_dx_aiuc(d->pbc, d->xtest, d->xref[i], dx);
 631                 r2 = norm2(dx);
 632                 if (r2 <= d->cutoff2)
 633                 {
 634                     if (action(d, i, r2))
 635                     {
 636                         d->prevnbi = nbi;
 637                         d->prevcai = cai;
 638                         d->previ   = i;
 639                         return true;
 640                     }
 641                 }
 642             }
 643             d->exclind = 0;
 644             cai        = 0;
 645         }
 646     }
 647     else
 648     {
 649         i = d->previ + 1;
 650         for (; i < d->nref; ++i)
 651         {
 652             if (is_excluded(d, i))
 653             {
 654                 continue;
 655             }
 656             if (d->pbc)
 657             {
 658                 pbc_dx(d->pbc, d->xtest, d->xref[i], dx);
 659             }
 660             else
 661             {
 662                 rvec_sub(d->xtest, d->xref[i], dx);
 663             }
 664             r2 = norm2(dx);
 665             if (r2 <= d->cutoff2)
 666             {
 667                 if (action(d, i, r2))
 668                 {
 669                     d->previ = i;
 670                     return true;
 671                 }
 672             }
 673         }
 674     }
 675     return false;
 676 }
 677
 678 /*! \brief
 679  * Helper function to use with grid_search() to find the next neighbor.
 680  *
 681  * Simply breaks the loop on the first found neighbor.
 682  */
 683 static bool
 684 within_action(gmx_ana_nbsearch_t *d, int i, real r2)
 685 {
 686     return true;
 687 }
 688
 689 /*! \brief
 690  * Helper function to use with grid_search() to find the minimum distance.
 691  */
 692 static bool
 693 mindist_action(gmx_ana_nbsearch_t *d, int i, real r2)
 694 {
 695     d->cutoff2 = r2;
 696     d->previ   = i;
 697     return false;
 698 }
 699
 700 /*!
 701  * \param[in] d   Neighborhood search data structure.
 702  * \param[in] x   Test position.
 703  * \returns   true if \p x is within the cutoff of any reference position,
 704  *   false otherwise.
 705  */
 706 bool
 707 gmx_ana_nbsearch_is_within(gmx_ana_nbsearch_t *d, const rvec x)
 708 {
 709     grid_search_start(d, x);
 710     return grid_search(d, &within_action);
 711 }
 712
 713 /*!
 714  * \param[in] d   Neighborhood search data structure.
 715  * \param[in] p   Test positions.
 716  * \param[in] i   Use the i'th position in \p p for testing.
 717  * \returns   true if the test position is within the cutoff of any reference
 718  *   position, false otherwise.
 719  */
 720 bool
 721 gmx_ana_nbsearch_pos_is_within(gmx_ana_nbsearch_t *d, const gmx_ana_pos_t *p, int i)
 722 {
 723     return gmx_ana_nbsearch_is_within(d, p->x[i]);
 724 }
 725
 726 /*!
 727  * \param[in] d   Neighborhood search data structure.
 728  * \param[in] x   Test position.
 729  * \returns   The distance to the nearest reference position, or the cutoff
 730  *   value if there are no reference positions within the cutoff.
 731  */
 732 real
 733 gmx_ana_nbsearch_mindist(gmx_ana_nbsearch_t *d, const rvec x)
 734 {
 735     real mind;
 736
 737     grid_search_start(d, x);
 738     grid_search(d, &mindist_action);
 739     mind       = sqrt(d->cutoff2);
 740     d->cutoff2 = sqr(d->cutoff);
 741     return mind;
 742 }
 743
 744 /*!
 745  * \param[in] d   Neighborhood search data structure.
 746  * \param[in] p   Test positions.
 747  * \param[in] i   Use the i'th position in \p p for testing.
 748  * \returns   The distance to the nearest reference position, or the cutoff
 749  *   value if there are no reference positions within the cutoff.
 750  */
 751 real
 752 gmx_ana_nbsearch_pos_mindist(gmx_ana_nbsearch_t *d, const gmx_ana_pos_t *p, int i)
 753 {
 754     return gmx_ana_nbsearch_mindist(d, p->x[i]);
 755 }
 756
 757 /*!
 758  * \param[in]  d   Neighborhood search data structure.
 759  * \param[in]  x   Test positions.
 760  * \param[out] jp  Index of the reference position in the first pair.
 761  * \returns    true if there are positions within the cutoff.
 762  */
 763 bool
 764 gmx_ana_nbsearch_first_within(gmx_ana_nbsearch_t *d, const rvec x, int *jp)
 765 {
 766     grid_search_start(d, x);
 767     return gmx_ana_nbsearch_next_within(d, jp);
 768 }
 769
 770 /*!
 771  * \param[in]  d   Neighborhood search data structure.
 772  * \param[in]  p   Test positions.
 773  * \param[in]  i   Use the i'th position in \p p.
 774  * \param[out] jp  Index of the reference position in the first pair.
 775  * \returns    true if there are positions within the cutoff.
 776  */
 777 bool
 778 gmx_ana_nbsearch_pos_first_within(gmx_ana_nbsearch_t *d, const gmx_ana_pos_t *p,
 779                                   int i, int *jp)
 780 {
 781     return gmx_ana_nbsearch_first_within(d, p->x[i], jp);
 782 }
 783
 784 /*!
 785  * \param[in]  d   Neighborhood search data structure.
 786  * \param[out] jp  Index of the test position in the next pair.
 787  * \returns    true if there are positions within the cutoff.
 788  */
 789 bool
 790 gmx_ana_nbsearch_next_within(gmx_ana_nbsearch_t *d, int *jp)
 791 {
 792     if (grid_search(d, &within_action))
 793     {
 794         *jp = d->previ;
 795         return true;
 796     }
 797     *jp = -1;
 798     return false;
 799 }
 800
 801 namespace gmx
 802 {
 803
 804 namespace internal
 805 {
 806
 807 /********************************************************************
 808  * AnalysisNeighborhoodPairSearchImpl
 809  */
 810
 811 class AnalysisNeighborhoodPairSearchImpl
 812 {
 813     public:
 814         AnalysisNeighborhoodPairSearchImpl()
 815             : nb_(NULL), testIndex_(0)
 816         {
 817         }
 818
 819         gmx_ana_nbsearch_t                 *nb_;
 820         int                                 testIndex_;
 821 };
 822
 823 /********************************************************************
 824  * AnalysisNeighborhoodSearchImpl
 825  */
 826
 827 class AnalysisNeighborhoodSearchImpl
 828 {
 829     public:
 830         typedef AnalysisNeighborhoodPairSearch::ImplPointer
 831             PairSearchImplPointer;
 832
 833         AnalysisNeighborhoodSearchImpl()
 834             : nb_(NULL), pairSearch_(new AnalysisNeighborhoodPairSearchImpl)
 835         {
 836         }
 837         ~AnalysisNeighborhoodSearchImpl()
 838         {
 839             GMX_RELEASE_ASSERT(pairSearch_.unique(),
 840                                "Dangling AnalysisNeighborhoodPairSearch reference");
 841             if (nb_ != NULL)
 842             {
 843                 gmx_ana_nbsearch_free(nb_);
 844             }
 845         }
 846
 847         gmx_ana_nbsearch_t     *nb_;
 848         PairSearchImplPointer   pairSearch_;
 849 };
 850
 851 }   // namespace internal
 852
 853 /********************************************************************
 854  * AnalysisNeighborhood::Impl
 855  */
 856
 857 class AnalysisNeighborhood::Impl
 858 {
 859     public:
 860         typedef AnalysisNeighborhoodSearch::ImplPointer SearchImplPointer;
 861
 862         Impl()
 863             : search_(new internal::AnalysisNeighborhoodSearchImpl), bTryGrid_(false)
 864         {
 865         }
 866         ~Impl()
 867         {
 868             GMX_RELEASE_ASSERT(search_.unique(),
 869                                "Dangling AnalysisNeighborhoodSearch reference");
 870         }
 871
 872         SearchImplPointer       search_;
 873         bool                    bTryGrid_;
 874 };
 875
 876 /********************************************************************
 877  * AnalysisNeighborhood
 878  */
 879
 880 AnalysisNeighborhood::AnalysisNeighborhood()
 881     : impl_(new Impl)
 882 {
 883 }
 884
 885 AnalysisNeighborhood::~AnalysisNeighborhood()
 886 {
 887 }
 888
 889 void AnalysisNeighborhood::setCutoff(real cutoff)
 890 {
 891     GMX_RELEASE_ASSERT(impl_->search_->nb_ == NULL,
 892                        "Multiple calls to setCutoff() not currently supported");
 893     impl_->search_->nb_              = gmx_ana_nbsearch_create(cutoff, 0);
 894     impl_->search_->pairSearch_->nb_ = impl_->search_->nb_;
 895     impl_->bTryGrid_                 = impl_->search_->nb_->bTryGrid;
 896 }
 897
 898 void AnalysisNeighborhood::setMode(SearchMode mode)
 899 {
 900     GMX_RELEASE_ASSERT(impl_->search_->nb_ != NULL,
 901                        "setParameters() not called");
 902     // TODO: Actually implement eSearchMode_Grid.
 903     switch (mode)
 904     {
 905         case eSearchMode_Automatic:
 906             impl_->search_->nb_->bTryGrid = impl_->bTryGrid_;
 907             break;
 908         case eSearchMode_Simple:
 909             impl_->search_->nb_->bTryGrid = false;
 910             break;
 911         case eSearchMode_Grid:
 912             impl_->search_->nb_->bTryGrid = impl_->bTryGrid_;
 913             break;
 914     }
 915 }
 916
 917 AnalysisNeighborhood::SearchMode AnalysisNeighborhood::mode() const
 918 {
 919     GMX_RELEASE_ASSERT(impl_->search_->nb_ != NULL,
 920                        "setParameters() not called");
 921     if (!impl_->search_->nb_->bTryGrid)
 922     {
 923         return eSearchMode_Simple;
 924     }
 925     return eSearchMode_Automatic;
 926 }
 927
 928 AnalysisNeighborhoodSearch
 929 AnalysisNeighborhood::initSearch(const t_pbc *pbc, int n, const rvec x[])
 930 {
 931     GMX_RELEASE_ASSERT(impl_->search_->nb_ != NULL,
 932                        "setParameters() not called");
 933     GMX_RELEASE_ASSERT(impl_->search_.unique(),
 934                        "Multiple concurrent searches not currently supported");
 935     gmx_ana_nbsearch_init(impl_->search_->nb_, const_cast<t_pbc *>(pbc), n, x);
 936     return AnalysisNeighborhoodSearch(impl_->search_);
 937 }
 938
 939 AnalysisNeighborhoodSearch
 940 AnalysisNeighborhood::initSearch(const t_pbc *pbc, const gmx_ana_pos_t *p)
 941 {
 942     GMX_RELEASE_ASSERT(impl_->search_->nb_ != NULL,
 943                        "setParameters() not called");
 944     GMX_RELEASE_ASSERT(impl_->search_.unique(),
 945                        "Multiple concurrent searches not currently supported");
 946     gmx_ana_nbsearch_pos_init(impl_->search_->nb_, const_cast<t_pbc *>(pbc), p);
 947     return AnalysisNeighborhoodSearch(impl_->search_);
 948 }
 949
 950 /********************************************************************
 951  * AnalysisNeighborhoodSearch
 952  */
 953
 954 AnalysisNeighborhoodSearch::AnalysisNeighborhoodSearch()
 955 {
 956 }
 957
 958 AnalysisNeighborhoodSearch::AnalysisNeighborhoodSearch(const ImplPointer &impl)
 959     : impl_(impl)
 960 {
 961 }
 962
 963 void AnalysisNeighborhoodSearch::reset()
 964 {
 965     impl_.reset();
 966 }
 967
 968 AnalysisNeighborhood::SearchMode AnalysisNeighborhoodSearch::mode() const
 969 {
 970     GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
 971     return (impl_->nb_->bGrid
 972             ? AnalysisNeighborhood::eSearchMode_Grid
 973             : AnalysisNeighborhood::eSearchMode_Simple);
 974 }
 975
 976 bool AnalysisNeighborhoodSearch::isWithin(const rvec x) const
 977 {
 978     GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
 979     GMX_RELEASE_ASSERT(impl_->pairSearch_.unique(),
 980                        "Calls to other methods concurrently with a pair search "
 981                        "not currently supported");
 982     return gmx_ana_nbsearch_is_within(impl_->nb_, x);
 983 }
 984
 985 bool AnalysisNeighborhoodSearch::isWithin(const gmx_ana_pos_t *p, int i) const
 986 {
 987     GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
 988     GMX_RELEASE_ASSERT(impl_->pairSearch_.unique(),
 989                        "Calls to other methods concurrently with a pair search "
 990                        "not currently supported");
 991     return gmx_ana_nbsearch_pos_is_within(impl_->nb_, p, i);
 992 }
 993
 994 real AnalysisNeighborhoodSearch::minimumDistance(const rvec x) const
 995 {
 996     GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
 997     GMX_RELEASE_ASSERT(impl_->pairSearch_.unique(),
 998                        "Calls to other methods concurrently with a pair search "
 999                        "not currently supported");
1000     return gmx_ana_nbsearch_mindist(impl_->nb_, x);
1001 }
1002
1003 real AnalysisNeighborhoodSearch::minimumDistance(const gmx_ana_pos_t *p, int i) const
1004 {
1005     GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
1006     GMX_RELEASE_ASSERT(impl_->pairSearch_.unique(),
1007                        "Calls to other methods concurrently with a pair search "
1008                        "not currently supported");
1009     return gmx_ana_nbsearch_pos_mindist(impl_->nb_, p, i);
1010 }
1011
1012 AnalysisNeighborhoodPair
1013 AnalysisNeighborhoodSearch::nearestPoint(const rvec x) const
1014 {
1015     GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
1016     GMX_RELEASE_ASSERT(impl_->pairSearch_.unique(),
1017                        "Calls to other methods concurrently with a pair search "
1018                        "not currently supported");
1019     (void)gmx_ana_nbsearch_mindist(impl_->nb_, x);
1020     return AnalysisNeighborhoodPair(impl_->nb_->previ, 0);
1021 }
1022
1023 AnalysisNeighborhoodPair
1024 AnalysisNeighborhoodSearch::nearestPoint(const gmx_ana_pos_t *p, int i) const
1025 {
1026     GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
1027     GMX_RELEASE_ASSERT(impl_->pairSearch_.unique(),
1028                        "Calls to other methods concurrently with a pair search "
1029                        "not currently supported");
1030     (void)gmx_ana_nbsearch_pos_mindist(impl_->nb_, p, i);
1031     return AnalysisNeighborhoodPair(impl_->nb_->previ, i);
1032 }
1033
1034 AnalysisNeighborhoodPairSearch
1035 AnalysisNeighborhoodSearch::startPairSearch(const rvec x)
1036 {
1037     GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
1038     GMX_RELEASE_ASSERT(impl_->pairSearch_.unique(),
1039                        "Multiple concurrent searches not currently supported");
1040     grid_search_start(impl_->nb_, x);
1041     impl_->pairSearch_->testIndex_ = 0;
1042     return AnalysisNeighborhoodPairSearch(impl_->pairSearch_);
1043 }
1044
1045 AnalysisNeighborhoodPairSearch
1046 AnalysisNeighborhoodSearch::startPairSearch(const gmx_ana_pos_t *p, int i)
1047 {
1048     GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
1049     GMX_RELEASE_ASSERT(impl_->pairSearch_.unique(),
1050                        "Multiple concurrent searches not currently supported");
1051     grid_search_start(impl_->nb_, p->x[i]);
1052     impl_->pairSearch_->testIndex_ = i;
1053     return AnalysisNeighborhoodPairSearch(impl_->pairSearch_);
1054 }
1055
1056 /********************************************************************
1057  * AnalysisNeighborhoodPairSearch
1058  */
1059
1060 AnalysisNeighborhoodPairSearch::AnalysisNeighborhoodPairSearch(
1061         const ImplPointer &impl)
1062     : impl_(impl)
1063 {
1064 }
1065
1066 bool AnalysisNeighborhoodPairSearch::findNextPair(AnalysisNeighborhoodPair *pair)
1067 {
1068     if (grid_search(impl_->nb_, &within_action))
1069     {
1070         *pair = AnalysisNeighborhoodPair(impl_->nb_->previ, impl_->testIndex_);
1071         return true;
1072     }
1073     *pair = AnalysisNeighborhoodPair();
1074     return false;
1075 }
1076
1077 } // namespace gmx