src/gromacs/selection/nbsearch.cpp

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