src/gmxlib/trajana/nbsearch.c

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