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