src/gromacs/selection/nbsearch.cpp

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