3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
9 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
10 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
11 * Copyright (c) 2001-2009, The GROMACS development team,
12 * check out http://www.gromacs.org for more information.
14 * This program is free software; you can redistribute it and/or
15 * modify it under the terms of the GNU General Public License
16 * as published by the Free Software Foundation; either version 2
17 * of the License, or (at your option) any later version.
19 * If you want to redistribute modifications, please consider that
20 * scientific software is very special. Version control is crucial -
21 * bugs must be traceable. We will be happy to consider code for
22 * inclusion in the official distribution, but derived work must not
23 * be called official GROMACS. Details are found in the README & COPYING
24 * files - if they are missing, get the official version at www.gromacs.org.
26 * To help us fund GROMACS development, we humbly ask that you cite
27 * the papers on the package - you can find them in the top README file.
29 * For more info, check our website at http://www.gromacs.org
31 /*! \page nbsearch Neighborhood search routines
33 * Functions to find particles within a neighborhood of a set of particles
34 * are defined in nbsearch.h.
35 * The usage is simple: a data structure is allocated with
36 * gmx_ana_nbsearch_create(), and the box shape and reference positions for a
37 * frame are set using gmx_ana_nbsearch_init() or gmx_ana_nbsearch_pos_init().
38 * Searches can then be performed with gmx_ana_nbsearch_is_within() and
39 * gmx_ana_nbsearch_mindist(), or with versions that take the \c gmx_ana_pos_t
41 * When the data structure is no longer required, it can be freed with
42 * gmx_ana_nbsearch_free().
47 * The grid implementation could still be optimized in several different ways:
48 * - Triclinic grid cells are not the most efficient shape, but make PBC
50 * - Precalculating the required PBC shift for a pair of cells outside the
51 * inner loop. After this is done, it should be quite straightforward to
52 * move to rectangular cells.
53 * - Pruning grid cells from the search list if they are completely outside
54 * the sphere that is being considered.
55 * - A better heuristic could be added for falling back to simple loops for a
56 * small number of reference particles.
57 * - A better heuristic for selecting the grid size.
58 * - A multi-level grid implementation could be used to be able to use small
59 * grids for short cutoffs with very inhomogeneous particle distributions
60 * without a memory cost.
64 * Implements functions in nbsearch.h.
66 * \author Teemu Murtola <teemu.murtola@cbr.su.se>
67 * \ingroup module_selection
75 #include "gromacs/legacyheaders/smalloc.h"
76 #include "gromacs/legacyheaders/typedefs.h"
77 #include "gromacs/legacyheaders/pbc.h"
78 #include "gromacs/legacyheaders/vec.h"
80 #include "gromacs/selection/nbsearch.h"
81 #include "gromacs/selection/position.h"
84 * Data structure for neighborhood searches.
86 struct gmx_ana_nbsearch_t
90 /** The cutoff squared. */
92 /** Maximum number of reference points. */
95 /** Number of reference points for the current frame. */
97 /** Reference point positions. */
99 /** Reference position ids (NULL if not available). */
104 /** Number of excluded reference positions for current test particle. */
106 /** Exclusions for current test particle. */
109 /** Whether to try grid searching. */
111 /** Whether grid searching is actually used for the current positions. */
113 /** Array allocated for storing in-unit-cell reference positions. */
115 /** false if the box is rectangular. */
117 /** Box vectors of a single grid cell. */
119 /** The reciprocal cell vectors as columns; the inverse of \p cellbox. */
121 /** Number of cells along each dimension. */
123 /** Total number of cells. */
125 /** Number of reference positions in each cell. */
127 /** List of reference positions in each cell. */
129 /** Allocation counts for each \p catom[i]. */
131 /** Allocation count for the per-cell arrays. */
133 /** Number of neighboring cells to consider. */
135 /** Offsets of the neighboring cells to consider. */
137 /** Allocation count for \p gnboffs. */
140 /** Stores test position during a pair loop. */
142 /** Stores the previous returned position during a pair loop. */
144 /** Stores the current exclusion index during loops. */
146 /** Stores the test particle cell index during loops. */
148 /** Stores the current cell neighbor index during pair loops. */
150 /** Stores the index within the current cell during pair loops. */
155 * \param[in] cutoff Cutoff distance for the search
156 * (<=0 stands for no cutoff).
157 * \param[in] maxn Maximum number of reference particles.
158 * \returns Created neighborhood search data structure.
161 gmx_ana_nbsearch_create(real cutoff, int maxn)
163 gmx_ana_nbsearch_t *d;
169 cutoff = GMX_REAL_MAX;
173 d->cutoff2 = sqr(cutoff);
180 d->xref_alloc = NULL;
189 d->gnboffs_nalloc = 0;
195 * \param d Data structure to free.
197 * After the call, the pointer \p d is no longer valid.
200 gmx_ana_nbsearch_free(gmx_ana_nbsearch_t *d)
202 sfree(d->xref_alloc);
208 for (ci = 0; ci < d->ncells; ++ci)
214 sfree(d->catom_nalloc);
220 * Calculates offsets to neighboring grid cells that should be considered.
222 * \param[in,out] d Grid information.
223 * \param[in] pbc Information about the box.
226 grid_init_cell_nblist(gmx_ana_nbsearch_t *d, t_pbc *pbc)
228 int maxx, maxy, maxz;
232 /* Find the extent of the sphere in triclinic coordinates */
233 maxz = (int)(d->cutoff * d->recipcell[ZZ][ZZ]) + 1;
234 rvnorm = sqrt(sqr(d->recipcell[YY][YY]) + sqr(d->recipcell[ZZ][YY]));
235 maxy = (int)(d->cutoff * rvnorm) + 1;
236 rvnorm = sqrt(sqr(d->recipcell[XX][XX]) + sqr(d->recipcell[YY][XX])
237 + sqr(d->recipcell[ZZ][XX]));
238 maxx = (int)(d->cutoff * rvnorm) + 1;
240 /* Calculate the number of cells and reallocate if necessary */
241 d->ngridnb = (2 * maxx + 1) * (2 * maxy + 1) * (2 * maxz + 1);
242 if (d->gnboffs_nalloc < d->ngridnb)
244 d->gnboffs_nalloc = d->ngridnb;
245 srenew(d->gnboffs, d->gnboffs_nalloc);
248 /* Store the whole cube */
249 /* TODO: Prune off corners that are not needed */
251 for (x = -maxx; x <= maxx; ++x)
253 for (y = -maxy; y <= maxy; ++y)
255 for (z = -maxz; z <= maxz; ++z)
257 d->gnboffs[i][XX] = x;
258 d->gnboffs[i][YY] = y;
259 d->gnboffs[i][ZZ] = z;
267 * Determines a suitable grid size.
269 * \param[in,out] d Grid information.
270 * \param[in] pbc Information about the box.
271 * \returns false if grid search is not suitable.
274 grid_setup_cells(gmx_ana_nbsearch_t *d, t_pbc *pbc)
280 targetsize = cbrt(pbc->box[XX][XX] * pbc->box[YY][YY] * pbc->box[ZZ][ZZ]
283 targetsize = pow(pbc->box[XX][XX] * pbc->box[YY][YY] * pbc->box[ZZ][ZZ]
284 * 10 / d->nref, (real)(1./3.));
288 for (dd = 0; dd < DIM; ++dd)
290 d->ncelldim[dd] = (int)(pbc->box[dd][dd] / targetsize);
291 d->ncells *= d->ncelldim[dd];
292 if (d->ncelldim[dd] < 3)
297 /* Reallocate if necessary */
298 if (d->cells_nalloc < d->ncells)
302 srenew(d->ncatoms, d->ncells);
303 srenew(d->catom, d->ncells);
304 srenew(d->catom_nalloc, d->ncells);
305 for (i = d->cells_nalloc; i < d->ncells; ++i)
308 d->catom_nalloc[i] = 0;
310 d->cells_nalloc = d->ncells;
316 * Sets ua a search grid for a given box.
318 * \param[in,out] d Grid information.
319 * \param[in] pbc Information about the box.
320 * \returns false if grid search is not suitable.
323 grid_set_box(gmx_ana_nbsearch_t *d, t_pbc *pbc)
327 /* TODO: This check could be improved. */
328 if (0.5*pbc->max_cutoff2 < d->cutoff2)
333 if (!grid_setup_cells(d, pbc))
338 d->bTric = TRICLINIC(pbc->box);
341 for (dd = 0; dd < DIM; ++dd)
343 svmul(1.0 / d->ncelldim[dd], pbc->box[dd], d->cellbox[dd]);
345 m_inv_ur0(d->cellbox, d->recipcell);
349 for (dd = 0; dd < DIM; ++dd)
351 d->cellbox[dd][dd] = pbc->box[dd][dd] / d->ncelldim[dd];
352 d->recipcell[dd][dd] = 1 / d->cellbox[dd][dd];
355 grid_init_cell_nblist(d, pbc);
360 * Maps a point into a grid cell.
362 * \param[in] d Grid information.
363 * \param[in] x Point to map.
364 * \param[out] cell Indices of the grid cell in which \p x lies.
366 * \p x should be in the triclinic unit cell.
369 grid_map_onto(gmx_ana_nbsearch_t *d, const rvec x, ivec cell)
377 tmvmul_ur0(d->recipcell, x, tx);
378 for (dd = 0; dd < DIM; ++dd)
380 cell[dd] = (int)tx[dd];
385 for (dd = 0; dd < DIM; ++dd)
387 cell[dd] = (int)(x[dd] * d->recipcell[dd][dd]);
393 * Calculates linear index of a grid cell.
395 * \param[in] d Grid information.
396 * \param[in] cell Cell indices.
397 * \returns Linear index of \p cell.
400 grid_index(gmx_ana_nbsearch_t *d, const ivec cell)
402 return cell[XX] + cell[YY] * d->ncelldim[XX]
403 + cell[ZZ] * d->ncelldim[XX] * d->ncelldim[YY];
407 * Clears all grid cells.
409 * \param[in,out] d Grid information.
412 grid_clear_cells(gmx_ana_nbsearch_t *d)
416 for (ci = 0; ci < d->ncells; ++ci)
423 * Adds an index into a grid cell.
425 * \param[in,out] d Grid information.
426 * \param[in] cell Cell into which \p i should be added.
427 * \param[in] i Index to add.
430 grid_add_to_cell(gmx_ana_nbsearch_t *d, const ivec cell, int i)
432 int ci = grid_index(d, cell);
434 if (d->ncatoms[ci] == d->catom_nalloc[ci])
436 d->catom_nalloc[ci] += 10;
437 srenew(d->catom[ci], d->catom_nalloc[ci]);
439 d->catom[ci][d->ncatoms[ci]++] = i;
443 * \param[in,out] d Neighborhood search data structure.
444 * \param[in] pbc PBC information for the frame.
445 * \param[in] n Number of reference positions for the frame.
446 * \param[in] x \p n reference positions for the frame.
448 * Initializes the data structure \p d such that it can be used to search
449 * for the neighbors of \p x.
452 gmx_ana_nbsearch_init(gmx_ana_nbsearch_t *d, t_pbc *pbc, int n, const rvec x[])
460 else if (d->bTryGrid)
462 d->bGrid = grid_set_box(d, pbc);
470 snew(d->xref_alloc, d->maxnref);
472 d->xref = d->xref_alloc;
475 for (i = 0; i < n; ++i)
477 copy_rvec(x[i], d->xref[i]);
479 put_atoms_in_triclinic_unitcell(ecenterTRIC, pbc->box, n, d->xref);
480 for (i = 0; i < n; ++i)
484 grid_map_onto(d, d->xref[i], refcell);
485 grid_add_to_cell(d, refcell, i);
490 // Won't be modified in this case, but when a grid is used,
491 // xref _is_ modified, so it can't be const.
492 d->xref = const_cast<rvec *>(x);
498 * \param[in,out] d Neighborhood search data structure.
499 * \param[in] pbc PBC information for the frame.
500 * \param[in] p Reference positions for the frame.
502 * A convenience wrapper for gmx_ana_nbsearch_init().
505 gmx_ana_nbsearch_pos_init(gmx_ana_nbsearch_t *d, t_pbc *pbc, const gmx_ana_pos_t *p)
507 gmx_ana_nbsearch_init(d, pbc, p->nr, p->x);
508 d->refid = (p->nr < d->maxnref ? p->m.refid : NULL);
512 * \param[in,out] d Neighborhood search data structure.
513 * \param[in] nexcl Number of reference positions to exclude from next
515 * \param[in] excl Indices of reference positions to exclude.
517 * The set exclusions remain in effect until the next call of this function.
520 gmx_ana_nbsearch_set_excl(gmx_ana_nbsearch_t *d, int nexcl, int excl[])
528 * Helper function to check whether a reference point should be excluded.
531 is_excluded(gmx_ana_nbsearch_t *d, int j)
533 if (d->exclind < d->nexcl)
537 while (d->exclind < d->nexcl && d->refid[j] > d->excl[d->exclind])
541 if (d->exclind < d->nexcl && d->refid[j] == d->excl[d->exclind])
549 while (d->bGrid && d->exclind < d->nexcl && d->excl[d->exclind] < j)
553 if (d->excl[d->exclind] == j)
564 * Initializes a grid search to find reference positions neighboring \p x.
567 grid_search_start(gmx_ana_nbsearch_t *d, const rvec x)
569 copy_rvec(x, d->xtest);
572 put_atoms_in_triclinic_unitcell(ecenterTRIC, d->pbc->box, 1, &d->xtest);
573 grid_map_onto(d, d->xtest, d->testcell);
585 * Does a grid search.
588 grid_search(gmx_ana_nbsearch_t *d,
589 bool (*action)(gmx_ana_nbsearch_t *d, int i, real r2))
600 cai = d->prevcai + 1;
602 for ( ; nbi < d->ngridnb; ++nbi)
606 ivec_add(d->testcell, d->gnboffs[nbi], cell);
607 /* TODO: Support for 2D and screw PBC */
608 cell[XX] = (cell[XX] + d->ncelldim[XX]) % d->ncelldim[XX];
609 cell[YY] = (cell[YY] + d->ncelldim[YY]) % d->ncelldim[YY];
610 cell[ZZ] = (cell[ZZ] + d->ncelldim[ZZ]) % d->ncelldim[ZZ];
611 ci = grid_index(d, cell);
612 /* TODO: Calculate the required PBC shift outside the inner loop */
613 for ( ; cai < d->ncatoms[ci]; ++cai)
615 i = d->catom[ci][cai];
616 if (is_excluded(d, i))
620 pbc_dx_aiuc(d->pbc, d->xtest, d->xref[i], dx);
622 if (r2 <= d->cutoff2)
624 if (action(d, i, r2))
640 for ( ; i < d->nref; ++i)
642 if (is_excluded(d, i))
648 pbc_dx(d->pbc, d->xtest, d->xref[i], dx);
652 rvec_sub(d->xtest, d->xref[i], dx);
655 if (r2 <= d->cutoff2)
657 if (action(d, i, r2))
669 * Helper function to use with grid_search() to find the next neighbor.
671 * Simply breaks the loop on the first found neighbor.
674 within_action(gmx_ana_nbsearch_t *d, int i, real r2)
680 * Helper function to use with grid_search() to find the minimum distance.
683 mindist_action(gmx_ana_nbsearch_t *d, int i, real r2)
690 * \param[in] d Neighborhood search data structure.
691 * \param[in] x Test position.
692 * \returns true if \p x is within the cutoff of any reference position,
696 gmx_ana_nbsearch_is_within(gmx_ana_nbsearch_t *d, const rvec x)
698 grid_search_start(d, x);
699 return grid_search(d, &within_action);
703 * \param[in] d Neighborhood search data structure.
704 * \param[in] p Test positions.
705 * \param[in] i Use the i'th position in \p p for testing.
706 * \returns true if the test position is within the cutoff of any reference
707 * position, false otherwise.
710 gmx_ana_nbsearch_pos_is_within(gmx_ana_nbsearch_t *d, const gmx_ana_pos_t *p, int i)
712 return gmx_ana_nbsearch_is_within(d, p->x[i]);
716 * \param[in] d Neighborhood search data structure.
717 * \param[in] x Test position.
718 * \returns The distance to the nearest reference position, or the cutoff
719 * value if there are no reference positions within the cutoff.
722 gmx_ana_nbsearch_mindist(gmx_ana_nbsearch_t *d, const rvec x)
726 grid_search_start(d, x);
727 grid_search(d, &mindist_action);
728 mind = sqrt(d->cutoff2);
729 d->cutoff2 = sqr(d->cutoff);
734 * \param[in] d Neighborhood search data structure.
735 * \param[in] p Test positions.
736 * \param[in] i Use the i'th position in \p p for testing.
737 * \returns The distance to the nearest reference position, or the cutoff
738 * value if there are no reference positions within the cutoff.
741 gmx_ana_nbsearch_pos_mindist(gmx_ana_nbsearch_t *d, const gmx_ana_pos_t *p, int i)
743 return gmx_ana_nbsearch_mindist(d, p->x[i]);
747 * \param[in] d Neighborhood search data structure.
748 * \param[in] x Test positions.
749 * \param[out] jp Index of the reference position in the first pair.
750 * \returns true if there are positions within the cutoff.
753 gmx_ana_nbsearch_first_within(gmx_ana_nbsearch_t *d, const rvec x, int *jp)
755 grid_search_start(d, x);
756 return gmx_ana_nbsearch_next_within(d, jp);
760 * \param[in] d Neighborhood search data structure.
761 * \param[in] p Test positions.
762 * \param[in] i Use the i'th position in \p p.
763 * \param[out] jp Index of the reference position in the first pair.
764 * \returns true if there are positions within the cutoff.
767 gmx_ana_nbsearch_pos_first_within(gmx_ana_nbsearch_t *d, const gmx_ana_pos_t *p,
770 return gmx_ana_nbsearch_first_within(d, p->x[i], jp);
774 * \param[in] d Neighborhood search data structure.
775 * \param[out] jp Index of the test position in the next pair.
776 * \returns true if there are positions within the cutoff.
779 gmx_ana_nbsearch_next_within(gmx_ana_nbsearch_t *d, int *jp)
781 if (grid_search(d, &within_action))