src/gmxlib/trajana/nbsearch.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2009, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012,2013, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 /*! \page nbsearch Neighborhood search routines
  39  *
  40  * Functions to find particles within a neighborhood of a set of particles
  41  * are defined in nbsearch.h.
  42  * The usage is simple: a data structure is allocated with
  43  * gmx_ana_nbsearch_create(), and the box shape and reference positions for a
  44  * frame are set using gmx_ana_nbsearch_init() or gmx_ana_nbsearch_pos_init().
  45  * Searches can then be performed with gmx_ana_nbsearch_is_within() and
  46  * gmx_ana_nbsearch_mindist(), or with versions that take the \c gmx_ana_pos_t
  47  * data structure.
  48  * When the data structure is no longer required, it can be freed with
  49  * gmx_ana_nbsearch_free().
  50  *
  51  * \internal
  52  *
  53  * \todo
  54  * The grid implementation could still be optimized in several different ways:
  55  *   - Triclinic grid cells are not the most efficient shape, but make PBC
  56  *     handling easier.
  57  *   - Precalculating the required PBC shift for a pair of cells outside the
  58  *     inner loop. After this is done, it should be quite straightforward to
  59  *     move to rectangular cells.
  60  *   - Pruning grid cells from the search list if they are completely outside
  61  *     the sphere that is being considered.
  62  *   - A better heuristic could be added for falling back to simple loops for a
  63  *     small number of reference particles.
  64  *   - A better heuristic for selecting the grid size.
  65  *   - A multi-level grid implementation could be used to be able to use small
  66  *     grids for short cutoffs with very inhomogeneous particle distributions
  67  *     without a memory cost.
  68  */
  69 /*! \internal \file
  70  * \brief Implementation of functions in nbsearch.h.
  71  */
  72 #ifdef HAVE_CONFIG_H
  73 #include <config.h>
  74 #endif
  75
  76 #include <math.h>
  77
  78 #include <smalloc.h>
  79 #include <typedefs.h>
  80 #include <pbc.h>
  81 #include <vec.h>
  82
  83 #include <nbsearch.h>
  84 #include <position.h>
  85
  86 /*! \internal \brief
  87  * Data structure for neighborhood searches.
  88  */
  89 struct gmx_ana_nbsearch_t
  90 {
  91     /** The cutoff. */
  92     real           cutoff;
  93     /** The cutoff squared. */
  94     real           cutoff2;
  95     /** Maximum number of reference points. */
  96     int            maxnref;
  97
  98     /** Number of reference points for the current frame. */
  99     int            nref;
 100     /** Reference point positions. */
 101     rvec          *xref;
 102     /** Reference position ids (NULL if not available). */
 103     int           *refid;
 104     /** PBC data. */
 105     t_pbc         *pbc;
 106
 107     /** Number of excluded reference positions for current test particle. */
 108     int            nexcl;
 109     /** Exclusions for current test particle. */
 110     int           *excl;
 111
 112     /** Whether to try grid searching. */
 113     gmx_bool           bTryGrid;
 114     /** Whether grid searching is actually used for the current positions. */
 115     gmx_bool           bGrid;
 116     /** Array allocated for storing in-unit-cell reference positions. */
 117     rvec          *xref_alloc;
 118     /** FALSE if the box is rectangular. */
 119     gmx_bool           bTric;
 120     /** Box vectors of a single grid cell. */
 121     matrix         cellbox;
 122     /** The reciprocal cell vectors as columns; the inverse of \p cellbox. */
 123     matrix         recipcell;
 124     /** Number of cells along each dimension. */
 125     ivec           ncelldim;
 126     /** Total number of cells. */
 127     int            ncells;
 128     /** Number of reference positions in each cell. */
 129     int           *ncatoms;
 130     /** List of reference positions in each cell. */
 131     atom_id      **catom;
 132     /** Allocation counts for each \p catom[i]. */
 133     int           *catom_nalloc;
 134     /** Allocation count for the per-cell arrays. */
 135     int            cells_nalloc;
 136     /** Number of neighboring cells to consider. */
 137     int            ngridnb;
 138     /** Offsets of the neighboring cells to consider. */
 139     ivec          *gnboffs;
 140     /** Allocation count for \p gnboffs. */
 141     int            gnboffs_nalloc;
 142
 143     /** Stores test position during a pair loop. */
 144     rvec           xtest;
 145     /** Stores the previous returned position during a pair loop. */
 146     int            previ;
 147     /** Stores the current exclusion index during loops. */
 148     int            exclind;
 149     /** Stores the test particle cell index during loops. */
 150     ivec           testcell;
 151     /** Stores the current cell neighbor index during pair loops. */
 152     int            prevnbi;
 153     /** Stores the index within the current cell during pair loops. */
 154     int            prevcai;
 155 };
 156
 157 /*!
 158  * \param[out] data   Neighborhood search data structure pointer to initialize.
 159  * \param[in]  cutoff Cutoff distance for the search
 160  *   (<=0 stands for no cutoff).
 161  * \param[in]  maxn   Maximum number of reference particles.
 162  * \returns    0 on success.
 163  */
 164 int
 165 gmx_ana_nbsearch_create(gmx_ana_nbsearch_t **data, real cutoff, int maxn)
 166 {
 167     gmx_ana_nbsearch_t *d;
 168
 169     snew(d, 1);
 170     d->bTryGrid = TRUE;
 171     if (cutoff <= 0)
 172     {
 173         cutoff = GMX_REAL_MAX;
 174         d->bTryGrid = FALSE;
 175     }
 176     d->cutoff = cutoff;
 177     d->cutoff2 = sqr(cutoff);
 178     d->maxnref = maxn;
 179
 180     d->xref = NULL;
 181     d->nexcl = 0;
 182     d->exclind = 0;
 183
 184     d->xref_alloc = NULL;
 185     d->ncells = 0;
 186     d->ncatoms = NULL;
 187     d->catom = NULL;
 188     d->catom_nalloc = 0;
 189     d->cells_nalloc = 0;
 190
 191     d->ngridnb = 0;
 192     d->gnboffs = NULL;
 193     d->gnboffs_nalloc = 0;
 194
 195     *data = d;
 196     return 0;
 197 }
 198
 199 /*!
 200  * \param     d Data structure to free.
 201  *
 202  * After the call, the pointer \p d is no longer valid.
 203  */
 204 void
 205 gmx_ana_nbsearch_free(gmx_ana_nbsearch_t *d)
 206 {
 207     sfree(d->xref_alloc);
 208     sfree(d->ncatoms);
 209     if (d->catom)
 210     {
 211         int ci;
 212
 213         for (ci = 0; ci < d->ncells; ++ci)
 214         {
 215             sfree(d->catom[ci]);
 216         }
 217         sfree(d->catom);
 218     }
 219     sfree(d->catom_nalloc);
 220     sfree(d->gnboffs);
 221     sfree(d);
 222 }
 223
 224 /*! \brief
 225  * Calculates offsets to neighboring grid cells that should be considered.
 226  *
 227  * \param[in,out] d    Grid information.
 228  * \param[in]     pbc  Information about the box.
 229  */
 230 static void
 231 grid_init_cell_nblist(gmx_ana_nbsearch_t *d, t_pbc *pbc)
 232 {
 233     int   maxx, maxy, maxz;
 234     int   x, y, z, i;
 235     real  rvnorm;
 236
 237     /* Find the extent of the sphere in triclinic coordinates */
 238     maxz = (int)(d->cutoff * d->recipcell[ZZ][ZZ]) + 1;
 239     rvnorm = sqrt(sqr(d->recipcell[YY][YY]) + sqr(d->recipcell[ZZ][YY]));
 240     maxy = (int)(d->cutoff * rvnorm) + 1;
 241     rvnorm = sqrt(sqr(d->recipcell[XX][XX]) + sqr(d->recipcell[YY][XX])
 242                   + sqr(d->recipcell[ZZ][XX]));
 243     maxx = (int)(d->cutoff * rvnorm) + 1;
 244
 245     /* Calculate the number of cells and reallocate if necessary */
 246     d->ngridnb = (2 * maxx + 1) * (2 * maxy + 1) * (2 * maxz + 1);
 247     if (d->gnboffs_nalloc < d->ngridnb)
 248     {
 249         d->gnboffs_nalloc = d->ngridnb;
 250         srenew(d->gnboffs, d->gnboffs_nalloc);
 251     }
 252
 253     /* Store the whole cube */
 254     /* TODO: Prune off corners that are not needed */
 255     i = 0;
 256     for (x = -maxx; x <= maxx; ++x)
 257     {
 258         for (y = -maxy; y <= maxy; ++y)
 259         {
 260             for (z = -maxz; z <= maxz; ++z)
 261             {
 262                 d->gnboffs[i][XX] = x;
 263                 d->gnboffs[i][YY] = y;
 264                 d->gnboffs[i][ZZ] = z;
 265                 ++i;
 266             }
 267         }
 268     }
 269 }
 270
 271 /*! \brief
 272  * Determines a suitable grid size.
 273  *
 274  * \param[in,out] d    Grid information.
 275  * \param[in]     pbc  Information about the box.
 276  * \returns  FALSE if grid search is not suitable.
 277  */
 278 static gmx_bool
 279 grid_setup_cells(gmx_ana_nbsearch_t *d, t_pbc *pbc)
 280 {
 281     real targetsize;
 282     int  dd;
 283
 284 #ifdef HAVE_CBRT
 285     targetsize = cbrt(pbc->box[XX][XX] * pbc->box[YY][YY] * pbc->box[ZZ][ZZ]
 286                       * 10 / d->nref);
 287 #else
 288     targetsize = pow(pbc->box[XX][XX] * pbc->box[YY][YY] * pbc->box[ZZ][ZZ]
 289                       * 10 / d->nref, 1./3.);
 290 #endif
 291
 292     d->ncells = 1;
 293     for (dd = 0; dd < DIM; ++dd)
 294     {
 295         d->ncelldim[dd] = (int)(pbc->box[dd][dd] / targetsize);
 296         d->ncells *= d->ncelldim[dd];
 297         if (d->ncelldim[dd] < 3)
 298         {
 299             return FALSE;
 300         }
 301     }
 302     /* Reallocate if necessary */
 303     if (d->cells_nalloc < d->ncells)
 304     {
 305         int  i;
 306
 307         srenew(d->ncatoms, d->ncells);
 308         srenew(d->catom, d->ncells);
 309         srenew(d->catom_nalloc, d->ncells);
 310         for (i = d->cells_nalloc; i < d->ncells; ++i)
 311         {
 312             d->catom[i] = NULL;
 313             d->catom_nalloc[i] = 0;
 314         }
 315         d->cells_nalloc = d->ncells;
 316     }
 317     return TRUE;
 318 }
 319
 320 /*! \brief
 321  * Sets ua a search grid for a given box.
 322  *
 323  * \param[in,out] d    Grid information.
 324  * \param[in]     pbc  Information about the box.
 325  * \returns  FALSE if grid search is not suitable.
 326  */
 327 static gmx_bool
 328 grid_set_box(gmx_ana_nbsearch_t *d, t_pbc *pbc)
 329 {
 330     int dd;
 331
 332     /* TODO: This check could be improved. */
 333     if (0.5*pbc->max_cutoff2 < d->cutoff2)
 334     {
 335         return FALSE;
 336     }
 337
 338     if (!grid_setup_cells(d, pbc))
 339     {
 340         return FALSE;
 341     }
 342
 343     d->bTric = TRICLINIC(pbc->box);
 344     if (d->bTric)
 345     {
 346         for (dd = 0; dd < DIM; ++dd)
 347         {
 348             svmul(1.0 / d->ncelldim[dd], pbc->box[dd], d->cellbox[dd]);
 349         }
 350         m_inv_ur0(d->cellbox, d->recipcell);
 351     }
 352     else
 353     {
 354         for (dd = 0; dd < DIM; ++dd)
 355         {
 356             d->cellbox[dd][dd] = pbc->box[dd][dd] / d->ncelldim[dd];
 357             d->recipcell[dd][dd] = 1 / d->cellbox[dd][dd];
 358         }
 359     }
 360     grid_init_cell_nblist(d, pbc);
 361     return TRUE;
 362 }
 363
 364 /*! \brief
 365  * Maps a point into a grid cell.
 366  *
 367  * \param[in]  d    Grid information.
 368  * \param[in]  x    Point to map.
 369  * \param[out] cell Indices of the grid cell in which \p x lies.
 370  *
 371  * \p x should be in the triclinic unit cell.
 372  */
 373 static void
 374 grid_map_onto(gmx_ana_nbsearch_t *d, const rvec x, ivec cell)
 375 {
 376     int dd;
 377
 378     if (d->bTric)
 379     {
 380         rvec tx;
 381
 382         tmvmul_ur0(d->recipcell, x, tx);
 383         for (dd = 0; dd < DIM; ++dd)
 384         {
 385             cell[dd] = (int)tx[dd];
 386         }
 387     }
 388     else
 389     {
 390         for (dd = 0; dd < DIM; ++dd)
 391         {
 392             cell[dd] = (int)(x[dd] * d->recipcell[dd][dd]);
 393         }
 394     }
 395 }
 396
 397 /*! \brief
 398  * Calculates linear index of a grid cell.
 399  *
 400  * \param[in]  d    Grid information.
 401  * \param[in]  cell Cell indices.
 402  * \returns    Linear index of \p cell.
 403  */
 404 static int
 405 grid_index(gmx_ana_nbsearch_t *d, const ivec cell)
 406 {
 407     return cell[XX] + cell[YY] * d->ncelldim[XX]
 408         + cell[ZZ] * d->ncelldim[XX] * d->ncelldim[YY];
 409 }
 410
 411 /*! \brief
 412  * Clears all grid cells.
 413  *
 414  * \param[in,out] d    Grid information.
 415  */
 416 static void
 417 grid_clear_cells(gmx_ana_nbsearch_t *d)
 418 {
 419     int  ci;
 420
 421     for (ci = 0; ci < d->ncells; ++ci)
 422     {
 423         d->ncatoms[ci] = 0;
 424     }
 425 }
 426
 427 /*! \brief
 428  * Adds an index into a grid cell.
 429  *
 430  * \param[in,out] d    Grid information.
 431  * \param[in]     cell Cell into which \p i should be added.
 432  * \param[in]     i    Index to add.
 433  */
 434 static void
 435 grid_add_to_cell(gmx_ana_nbsearch_t *d, const ivec cell, int i)
 436 {
 437     int ci = grid_index(d, cell);
 438
 439     if (d->ncatoms[ci] == d->catom_nalloc[ci])
 440     {
 441         d->catom_nalloc[ci] += 10;
 442         srenew(d->catom[ci], d->catom_nalloc[ci]);
 443     }
 444     d->catom[ci][d->ncatoms[ci]++] = i;
 445 }
 446
 447 /*!
 448  * \param[in,out] d   Neighborhood search data structure.
 449  * \param[in]     pbc PBC information for the frame.
 450  * \param[in]     n   Number of reference positions for the frame.
 451  * \param[in]     x   \p n reference positions for the frame.
 452  * \returns       0 on success.
 453  *
 454  * Initializes the data structure \p d such that it can be used to search
 455  * for the neighbors of \p x.
 456  */
 457 int
 458 gmx_ana_nbsearch_init(gmx_ana_nbsearch_t *d, t_pbc *pbc, int n, rvec x[])
 459 {
 460     d->pbc  = pbc;
 461     d->nref = n;
 462     if (!pbc)
 463     {
 464         d->bGrid = FALSE;
 465     }
 466     else if (d->bTryGrid)
 467     {
 468         d->bGrid = grid_set_box(d, pbc);
 469     }
 470     if (d->bGrid)
 471     {
 472         int  i;
 473
 474         if (!d->xref_alloc)
 475         {
 476             snew(d->xref_alloc, d->maxnref);
 477         }
 478         d->xref = d->xref_alloc;
 479         grid_clear_cells(d);
 480
 481         for (i = 0; i < n; ++i)
 482         {
 483             copy_rvec(x[i], d->xref[i]);
 484         }
 485         put_atoms_in_triclinic_unitcell(ecenterTRIC, pbc->box, n, d->xref);
 486         for (i = 0; i < n; ++i)
 487         {
 488             ivec refcell;
 489
 490             grid_map_onto(d, d->xref[i], refcell);
 491             grid_add_to_cell(d, refcell, i);
 492         }
 493     }
 494     else
 495     {
 496         d->xref = x;
 497     }
 498     d->refid = NULL;
 499     return 0;
 500 }
 501
 502 /*!
 503  * \param[in,out] d   Neighborhood search data structure.
 504  * \param[in]     pbc PBC information for the frame.
 505  * \param[in]     p   Reference positions for the frame.
 506  * \returns       0 on success.
 507  *
 508  * A convenience wrapper for gmx_ana_nbsearch_init().
 509  */
 510 int
 511 gmx_ana_nbsearch_pos_init(gmx_ana_nbsearch_t *d, t_pbc *pbc, gmx_ana_pos_t *p)
 512 {
 513     int rc;
 514
 515     rc = gmx_ana_nbsearch_init(d, pbc, p->nr, p->x);
 516     d->refid = (p->nr < d->maxnref ? p->m.refid : NULL);
 517     return rc;
 518 }
 519
 520 /*!
 521  * \param[in,out] d     Neighborhood search data structure.
 522  * \param[in]     nexcl Number of reference positions to exclude from next
 523  *      search.
 524  * \param[in]     excl  Indices of reference positions to exclude.
 525  * \returns       0 on success.
 526  *
 527  * The set exclusions remain in effect until the next call of this function.
 528  */
 529 int
 530 gmx_ana_nbsearch_set_excl(gmx_ana_nbsearch_t *d, int nexcl, int excl[])
 531 {
 532
 533     d->nexcl = nexcl;
 534     d->excl = excl;
 535     return 0;
 536 }
 537
 538 /*! \brief
 539  * Helper function to check whether a reference point should be excluded.
 540  */
 541 static gmx_bool
 542 is_excluded(gmx_ana_nbsearch_t *d, int j)
 543 {
 544     if (d->exclind < d->nexcl)
 545     {
 546         if (d->refid)
 547         {
 548             while (d->exclind < d->nexcl && d->refid[j] > d->excl[d->exclind])
 549             {
 550                 ++d->exclind;
 551             }
 552             if (d->exclind < d->nexcl && d->refid[j] == d->excl[d->exclind])
 553             {
 554                 ++d->exclind;
 555                 return TRUE;
 556             }
 557         }
 558         else
 559         {
 560             while (d->bGrid && d->exclind < d->nexcl && d->excl[d->exclind] < j)
 561             {
 562                 ++d->exclind;
 563             }
 564             if (d->excl[d->exclind] == j)
 565             {
 566                 ++d->exclind;
 567                 return TRUE;
 568             }
 569         }
 570     }
 571     return FALSE;
 572 }
 573
 574 /*! \brief
 575  * Initializes a grid search to find reference positions neighboring \p x.
 576  */
 577 static void
 578 grid_search_start(gmx_ana_nbsearch_t *d, rvec x)
 579 {
 580     copy_rvec(x, d->xtest);
 581     if (d->bGrid)
 582     {
 583         put_atoms_in_triclinic_unitcell(ecenterTRIC, d->pbc->box, 1, &d->xtest);
 584         grid_map_onto(d, d->xtest, d->testcell);
 585         d->prevnbi = 0;
 586         d->prevcai = -1;
 587     }
 588     else
 589     {
 590         d->previ = -1;
 591     }
 592     d->exclind = 0;
 593 }
 594
 595 /*! \brief
 596  * Does a grid search.
 597  */
 598 static gmx_bool
 599 grid_search(gmx_ana_nbsearch_t *d,
 600             gmx_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             /* TODO: Support for 2D and screw PBC */
 619             cell[XX] = (cell[XX] + d->ncelldim[XX]) % d->ncelldim[XX];
 620             cell[YY] = (cell[YY] + d->ncelldim[YY]) % d->ncelldim[YY];
 621             cell[ZZ] = (cell[ZZ] + d->ncelldim[ZZ]) % d->ncelldim[ZZ];
 622             ci = grid_index(d, cell);
 623             /* TODO: Calculate the required PBC shift outside the inner loop */
 624             for ( ; cai < d->ncatoms[ci]; ++cai)
 625             {
 626                 i = d->catom[ci][cai];
 627                 if (is_excluded(d, i))
 628                 {
 629                     continue;
 630                 }
 631                 pbc_dx_aiuc(d->pbc, d->xtest, d->xref[i], dx);
 632                 r2 = norm2(dx);
 633                 if (r2 <= d->cutoff2)
 634                 {
 635                     if (action(d, i, r2))
 636                     {
 637                         d->prevnbi = nbi;
 638                         d->prevcai = cai;
 639                         d->previ   = i;
 640                         return TRUE;
 641                     }
 642                 }
 643             }
 644             d->exclind = 0;
 645             cai = 0;
 646         }
 647     }
 648     else
 649     {
 650         i = d->previ + 1;
 651         for ( ; i < d->nref; ++i)
 652         {
 653             if (is_excluded(d, i))
 654             {
 655                 continue;
 656             }
 657             if (d->pbc)
 658             {
 659                 pbc_dx(d->pbc, d->xtest, d->xref[i], dx);
 660             }
 661             else
 662             {
 663                 rvec_sub(d->xtest, d->xref[i], dx);
 664             }
 665             r2 = norm2(dx);
 666             if (r2 <= d->cutoff2)
 667             {
 668                 if (action(d, i, r2))
 669                 {
 670                     d->previ = i;
 671                     return TRUE;
 672                 }
 673             }
 674         }
 675     }
 676     return FALSE;
 677 }
 678
 679 /*! \brief
 680  * Helper function to use with grid_search() to find the next neighbor.
 681  *
 682  * Simply breaks the loop on the first found neighbor.
 683  */
 684 static gmx_bool
 685 within_action(gmx_ana_nbsearch_t *d, int i, real r2)
 686 {
 687     return TRUE;
 688 }
 689
 690 /*! \brief
 691  * Helper function to use with grid_search() to find the minimum distance.
 692  */
 693 static gmx_bool
 694 mindist_action(gmx_ana_nbsearch_t *d, int i, real r2)
 695 {
 696     d->cutoff2 = r2;
 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 gmx_bool
 707 gmx_ana_nbsearch_is_within(gmx_ana_nbsearch_t *d, 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 gmx_bool
 721 gmx_ana_nbsearch_pos_is_within(gmx_ana_nbsearch_t *d, 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, 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, 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]  n   Number of test positions in \p x.
 760  * \param[in]  x   Test positions.
 761  * \param[out] jp  Index of the reference position in the first pair.
 762  * \returns    TRUE if there are positions within the cutoff.
 763  */
 764 gmx_bool
 765 gmx_ana_nbsearch_first_within(gmx_ana_nbsearch_t *d, rvec x, int *jp)
 766 {
 767     grid_search_start(d, x);
 768     return gmx_ana_nbsearch_next_within(d, jp);
 769 }
 770
 771 /*!
 772  * \param[in]  d   Neighborhood search data structure.
 773  * \param[in]  p   Test positions.
 774  * \param[in]  i   Use the i'th position in \p p.
 775  * \param[out] jp  Index of the reference position in the first pair.
 776  * \returns    TRUE if there are positions within the cutoff.
 777  */
 778 gmx_bool
 779 gmx_ana_nbsearch_pos_first_within(gmx_ana_nbsearch_t *d, gmx_ana_pos_t *p,
 780                                   int i, int *jp)
 781 {
 782     return gmx_ana_nbsearch_first_within(d, p->x[i], jp);
 783 }
 784
 785 /*!
 786  * \param[in]  d   Neighborhood search data structure.
 787  * \param[out] jp  Index of the test position in the next pair.
 788  * \returns    TRUE if there are positions within the cutoff.
 789  */
 790 gmx_bool
 791 gmx_ana_nbsearch_next_within(gmx_ana_nbsearch_t *d, int *jp)
 792 {
 793     if (grid_search(d, &within_action))
 794     {
 795         *jp = d->previ;
 796         return TRUE;
 797     }
 798     *jp = -1;
 799     return FALSE;
 800 }