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, 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 "gromacs/legacyheaders/smalloc.h"
  80 #include "gromacs/legacyheaders/typedefs.h"
  81 #include "gromacs/legacyheaders/pbc.h"
  82 #include "gromacs/legacyheaders/vec.h"
  83
  84 #include "gromacs/selection/nbsearch.h"
  85 #include "gromacs/selection/position.h"
  86
  87 /*! \internal \brief
  88  * Data structure for neighborhood searches.
  89  */
  90 struct gmx_ana_nbsearch_t
  91 {
  92     /** The cutoff. */
  93     real           cutoff;
  94     /** The cutoff squared. */
  95     real           cutoff2;
  96     /** Maximum number of reference points. */
  97     int            maxnref;
  98
  99     /** Number of reference points for the current frame. */
 100     int            nref;
 101     /** Reference point positions. */
 102     rvec          *xref;
 103     /** Reference position ids (NULL if not available). */
 104     int           *refid;
 105     /** PBC data. */
 106     t_pbc         *pbc;
 107
 108     /** Number of excluded reference positions for current test particle. */
 109     int            nexcl;
 110     /** Exclusions for current test particle. */
 111     int           *excl;
 112
 113     /** Whether to try grid searching. */
 114     bool           bTryGrid;
 115     /** Whether grid searching is actually used for the current positions. */
 116     bool           bGrid;
 117     /** Array allocated for storing in-unit-cell reference positions. */
 118     rvec          *xref_alloc;
 119     /** false if the box is rectangular. */
 120     bool           bTric;
 121     /** Box vectors of a single grid cell. */
 122     matrix         cellbox;
 123     /** The reciprocal cell vectors as columns; the inverse of \p cellbox. */
 124     matrix         recipcell;
 125     /** Number of cells along each dimension. */
 126     ivec           ncelldim;
 127     /** Total number of cells. */
 128     int            ncells;
 129     /** Number of reference positions in each cell. */
 130     int           *ncatoms;
 131     /** List of reference positions in each cell. */
 132     atom_id      **catom;
 133     /** Allocation counts for each \p catom[i]. */
 134     int           *catom_nalloc;
 135     /** Allocation count for the per-cell arrays. */
 136     int            cells_nalloc;
 137     /** Number of neighboring cells to consider. */
 138     int            ngridnb;
 139     /** Offsets of the neighboring cells to consider. */
 140     ivec          *gnboffs;
 141     /** Allocation count for \p gnboffs. */
 142     int            gnboffs_nalloc;
 143
 144     /** Stores test position during a pair loop. */
 145     rvec           xtest;
 146     /** Stores the previous returned position during a pair loop. */
 147     int            previ;
 148     /** Stores the current exclusion index during loops. */
 149     int            exclind;
 150     /** Stores the test particle cell index during loops. */
 151     ivec           testcell;
 152     /** Stores the current cell neighbor index during pair loops. */
 153     int            prevnbi;
 154     /** Stores the index within the current cell during pair loops. */
 155     int            prevcai;
 156 };
 157
 158 /*!
 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  Created neighborhood search data structure.
 163  */
 164 gmx_ana_nbsearch_t *
 165 gmx_ana_nbsearch_create(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     return d;
 196 }
 197
 198 /*!
 199  * \param     d Data structure to free.
 200  *
 201  * After the call, the pointer \p d is no longer valid.
 202  */
 203 void
 204 gmx_ana_nbsearch_free(gmx_ana_nbsearch_t *d)
 205 {
 206     sfree(d->xref_alloc);
 207     sfree(d->ncatoms);
 208     if (d->catom)
 209     {
 210         int ci;
 211
 212         for (ci = 0; ci < d->ncells; ++ci)
 213         {
 214             sfree(d->catom[ci]);
 215         }
 216         sfree(d->catom);
 217     }
 218     sfree(d->catom_nalloc);
 219     sfree(d->gnboffs);
 220     sfree(d);
 221 }
 222
 223 /*! \brief
 224  * Calculates offsets to neighboring grid cells that should be considered.
 225  *
 226  * \param[in,out] d    Grid information.
 227  * \param[in]     pbc  Information about the box.
 228  */
 229 static void
 230 grid_init_cell_nblist(gmx_ana_nbsearch_t *d, t_pbc *pbc)
 231 {
 232     int   maxx, maxy, maxz;
 233     int   x, y, z, i;
 234     real  rvnorm;
 235
 236     /* Find the extent of the sphere in triclinic coordinates */
 237     maxz   = (int)(d->cutoff * d->recipcell[ZZ][ZZ]) + 1;
 238     rvnorm = sqrt(sqr(d->recipcell[YY][YY]) + sqr(d->recipcell[ZZ][YY]));
 239     maxy   = (int)(d->cutoff * rvnorm) + 1;
 240     rvnorm = sqrt(sqr(d->recipcell[XX][XX]) + sqr(d->recipcell[YY][XX])
 241                   + sqr(d->recipcell[ZZ][XX]));
 242     maxx = (int)(d->cutoff * rvnorm) + 1;
 243
 244     /* Calculate the number of cells and reallocate if necessary */
 245     d->ngridnb = (2 * maxx + 1) * (2 * maxy + 1) * (2 * maxz + 1);
 246     if (d->gnboffs_nalloc < d->ngridnb)
 247     {
 248         d->gnboffs_nalloc = d->ngridnb;
 249         srenew(d->gnboffs, d->gnboffs_nalloc);
 250     }
 251
 252     /* Store the whole cube */
 253     /* TODO: Prune off corners that are not needed */
 254     i = 0;
 255     for (x = -maxx; x <= maxx; ++x)
 256     {
 257         for (y = -maxy; y <= maxy; ++y)
 258         {
 259             for (z = -maxz; z <= maxz; ++z)
 260             {
 261                 d->gnboffs[i][XX] = x;
 262                 d->gnboffs[i][YY] = y;
 263                 d->gnboffs[i][ZZ] = z;
 264                 ++i;
 265             }
 266         }
 267     }
 268 }
 269
 270 /*! \brief
 271  * Determines a suitable grid size.
 272  *
 273  * \param[in,out] d    Grid information.
 274  * \param[in]     pbc  Information about the box.
 275  * \returns  false if grid search is not suitable.
 276  */
 277 static bool
 278 grid_setup_cells(gmx_ana_nbsearch_t *d, t_pbc *pbc)
 279 {
 280     real targetsize;
 281     int  dd;
 282
 283 #ifdef HAVE_CBRT
 284     targetsize = cbrt(pbc->box[XX][XX] * pbc->box[YY][YY] * pbc->box[ZZ][ZZ]
 285                       * 10 / d->nref);
 286 #else
 287     targetsize = pow(pbc->box[XX][XX] * pbc->box[YY][YY] * pbc->box[ZZ][ZZ]
 288                      * 10 / d->nref, (real)(1./3.));
 289 #endif
 290
 291     d->ncells = 1;
 292     for (dd = 0; dd < DIM; ++dd)
 293     {
 294         d->ncelldim[dd] = (int)(pbc->box[dd][dd] / targetsize);
 295         d->ncells      *= d->ncelldim[dd];
 296         if (d->ncelldim[dd] < 3)
 297         {
 298             return false;
 299         }
 300     }
 301     /* Reallocate if necessary */
 302     if (d->cells_nalloc < d->ncells)
 303     {
 304         int  i;
 305
 306         srenew(d->ncatoms, d->ncells);
 307         srenew(d->catom, d->ncells);
 308         srenew(d->catom_nalloc, d->ncells);
 309         for (i = d->cells_nalloc; i < d->ncells; ++i)
 310         {
 311             d->catom[i]        = NULL;
 312             d->catom_nalloc[i] = 0;
 313         }
 314         d->cells_nalloc = d->ncells;
 315     }
 316     return true;
 317 }
 318
 319 /*! \brief
 320  * Sets ua a search grid for a given box.
 321  *
 322  * \param[in,out] d    Grid information.
 323  * \param[in]     pbc  Information about the box.
 324  * \returns  false if grid search is not suitable.
 325  */
 326 static bool
 327 grid_set_box(gmx_ana_nbsearch_t *d, t_pbc *pbc)
 328 {
 329     int dd;
 330
 331     /* TODO: This check could be improved. */
 332     if (0.5*pbc->max_cutoff2 < d->cutoff2)
 333     {
 334         return false;
 335     }
 336
 337     if (!grid_setup_cells(d, pbc))
 338     {
 339         return false;
 340     }
 341
 342     d->bTric = TRICLINIC(pbc->box);
 343     if (d->bTric)
 344     {
 345         for (dd = 0; dd < DIM; ++dd)
 346         {
 347             svmul(1.0 / d->ncelldim[dd], pbc->box[dd], d->cellbox[dd]);
 348         }
 349         m_inv_ur0(d->cellbox, d->recipcell);
 350     }
 351     else
 352     {
 353         for (dd = 0; dd < DIM; ++dd)
 354         {
 355             d->cellbox[dd][dd]   = pbc->box[dd][dd] / d->ncelldim[dd];
 356             d->recipcell[dd][dd] = 1 / d->cellbox[dd][dd];
 357         }
 358     }
 359     grid_init_cell_nblist(d, pbc);
 360     return true;
 361 }
 362
 363 /*! \brief
 364  * Maps a point into a grid cell.
 365  *
 366  * \param[in]  d    Grid information.
 367  * \param[in]  x    Point to map.
 368  * \param[out] cell Indices of the grid cell in which \p x lies.
 369  *
 370  * \p x should be in the triclinic unit cell.
 371  */
 372 static void
 373 grid_map_onto(gmx_ana_nbsearch_t *d, const rvec x, ivec cell)
 374 {
 375     int dd;
 376
 377     if (d->bTric)
 378     {
 379         rvec tx;
 380
 381         tmvmul_ur0(d->recipcell, x, tx);
 382         for (dd = 0; dd < DIM; ++dd)
 383         {
 384             cell[dd] = (int)tx[dd];
 385         }
 386     }
 387     else
 388     {
 389         for (dd = 0; dd < DIM; ++dd)
 390         {
 391             cell[dd] = (int)(x[dd] * d->recipcell[dd][dd]);
 392         }
 393     }
 394 }
 395
 396 /*! \brief
 397  * Calculates linear index of a grid cell.
 398  *
 399  * \param[in]  d    Grid information.
 400  * \param[in]  cell Cell indices.
 401  * \returns    Linear index of \p cell.
 402  */
 403 static int
 404 grid_index(gmx_ana_nbsearch_t *d, const ivec cell)
 405 {
 406     return cell[XX] + cell[YY] * d->ncelldim[XX]
 407            + cell[ZZ] * d->ncelldim[XX] * d->ncelldim[YY];
 408 }
 409
 410 /*! \brief
 411  * Clears all grid cells.
 412  *
 413  * \param[in,out] d    Grid information.
 414  */
 415 static void
 416 grid_clear_cells(gmx_ana_nbsearch_t *d)
 417 {
 418     int  ci;
 419
 420     for (ci = 0; ci < d->ncells; ++ci)
 421     {
 422         d->ncatoms[ci] = 0;
 423     }
 424 }
 425
 426 /*! \brief
 427  * Adds an index into a grid cell.
 428  *
 429  * \param[in,out] d    Grid information.
 430  * \param[in]     cell Cell into which \p i should be added.
 431  * \param[in]     i    Index to add.
 432  */
 433 static void
 434 grid_add_to_cell(gmx_ana_nbsearch_t *d, const ivec cell, int i)
 435 {
 436     int ci = grid_index(d, cell);
 437
 438     if (d->ncatoms[ci] == d->catom_nalloc[ci])
 439     {
 440         d->catom_nalloc[ci] += 10;
 441         srenew(d->catom[ci], d->catom_nalloc[ci]);
 442     }
 443     d->catom[ci][d->ncatoms[ci]++] = i;
 444 }
 445
 446 /*!
 447  * \param[in,out] d   Neighborhood search data structure.
 448  * \param[in]     pbc PBC information for the frame.
 449  * \param[in]     n   Number of reference positions for the frame.
 450  * \param[in]     x   \p n reference positions for the frame.
 451  *
 452  * Initializes the data structure \p d such that it can be used to search
 453  * for the neighbors of \p x.
 454  */
 455 void
 456 gmx_ana_nbsearch_init(gmx_ana_nbsearch_t *d, t_pbc *pbc, int n, const rvec x[])
 457 {
 458     d->pbc  = pbc;
 459     d->nref = n;
 460     if (!pbc)
 461     {
 462         d->bGrid = false;
 463     }
 464     else if (d->bTryGrid)
 465     {
 466         d->bGrid = grid_set_box(d, pbc);
 467     }
 468     if (d->bGrid)
 469     {
 470         int  i;
 471
 472         if (!d->xref_alloc)
 473         {
 474             snew(d->xref_alloc, d->maxnref);
 475         }
 476         d->xref = d->xref_alloc;
 477         grid_clear_cells(d);
 478
 479         for (i = 0; i < n; ++i)
 480         {
 481             copy_rvec(x[i], d->xref[i]);
 482         }
 483         put_atoms_in_triclinic_unitcell(ecenterTRIC, pbc->box, n, d->xref);
 484         for (i = 0; i < n; ++i)
 485         {
 486             ivec refcell;
 487
 488             grid_map_onto(d, d->xref[i], refcell);
 489             grid_add_to_cell(d, refcell, i);
 490         }
 491     }
 492     else
 493     {
 494         // Won't be modified in this case, but when a grid is used,
 495         // xref _is_ modified, so it can't be const.
 496         d->xref = const_cast<rvec *>(x);
 497     }
 498     d->refid = NULL;
 499 }
 500
 501 /*!
 502  * \param[in,out] d   Neighborhood search data structure.
 503  * \param[in]     pbc PBC information for the frame.
 504  * \param[in]     p   Reference positions for the frame.
 505  *
 506  * A convenience wrapper for gmx_ana_nbsearch_init().
 507  */
 508 void
 509 gmx_ana_nbsearch_pos_init(gmx_ana_nbsearch_t *d, t_pbc *pbc, const gmx_ana_pos_t *p)
 510 {
 511     gmx_ana_nbsearch_init(d, pbc, p->nr, p->x);
 512     d->refid = (p->nr < d->maxnref ? p->m.refid : NULL);
 513 }
 514
 515 /*!
 516  * \param[in,out] d     Neighborhood search data structure.
 517  * \param[in]     nexcl Number of reference positions to exclude from next
 518  *      search.
 519  * \param[in]     excl  Indices of reference positions to exclude.
 520  *
 521  * The set exclusions remain in effect until the next call of this function.
 522  */
 523 void
 524 gmx_ana_nbsearch_set_excl(gmx_ana_nbsearch_t *d, int nexcl, int excl[])
 525 {
 526
 527     d->nexcl = nexcl;
 528     d->excl  = excl;
 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, const 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, const 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, const 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, const 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, const 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]  x   Test positions.
 753  * \param[out] jp  Index of the reference position in the first pair.
 754  * \returns    true if there are positions within the cutoff.
 755  */
 756 bool
 757 gmx_ana_nbsearch_first_within(gmx_ana_nbsearch_t *d, const rvec x, int *jp)
 758 {
 759     grid_search_start(d, x);
 760     return gmx_ana_nbsearch_next_within(d, jp);
 761 }
 762
 763 /*!
 764  * \param[in]  d   Neighborhood search data structure.
 765  * \param[in]  p   Test positions.
 766  * \param[in]  i   Use the i'th position in \p p.
 767  * \param[out] jp  Index of the reference position in the first pair.
 768  * \returns    true if there are positions within the cutoff.
 769  */
 770 bool
 771 gmx_ana_nbsearch_pos_first_within(gmx_ana_nbsearch_t *d, const gmx_ana_pos_t *p,
 772                                   int i, int *jp)
 773 {
 774     return gmx_ana_nbsearch_first_within(d, p->x[i], jp);
 775 }
 776
 777 /*!
 778  * \param[in]  d   Neighborhood search data structure.
 779  * \param[out] jp  Index of the test position in the next pair.
 780  * \returns    true if there are positions within the cutoff.
 781  */
 782 bool
 783 gmx_ana_nbsearch_next_within(gmx_ana_nbsearch_t *d, int *jp)
 784 {
 785     if (grid_search(d, &within_action))
 786     {
 787         *jp = d->previ;
 788         return true;
 789     }
 790     *jp = -1;
 791     return false;
 792 }