index | modules | gitlab | page source

# Building neighbor-lists

A neighbor list is a collision detector for spheres: Given a number of spheres of different radius located at different points, it calculates the pairs of spheres that overlap.

ASE provides two implementations of neighbor lists. The newer linearly scaling function `neighbor_list()` and the older quadratically scaling class `PrimitiveNeighborList`. The latter will likely use the former as a backend in the future for linear scaling.

For flexibility, both implementations provide a “primitive” interface which accepts arrays as arguments rather than the more complex `Atoms` objects.

Both implementations can be used via the `NeighborList` class. It also provides easy access to the two implementations methods and functions. Constructing such an object can be done manually or with the `build_neighbor_list()` function.

Further functions provide access to some derived results like graph-analysis etc.:

> - `natural_cutoffs()`

> - `get_connectivity_matrix()`

> - `get_distance_matrix()`

> - `get_distance_indices()`

## API

class ase.neighborlist.NeighborList(cutoffs, skin=0.3, sorted=False, self_interaction=True, bothways=False, primitive=<class 'ase.neighborlist.PrimitiveNeighborList'>)

  Neighbor list object.

  cutoffs: list of float

    List of cutoff radii - one for each atom. If the spheres (defined by their cutoff radii) of two atoms overlap, they will be counted as neighbors. See `natural_cutoffs()` for an example on how to get such a list.

  skin: float

    If no atom has moved more than the skin-distance since the last call to the `update()` method, then the neighbor list can be reused. This will save some expensive rebuilds of the list, but extra neighbors outside the cutoff will be returned.

  self_interaction: bool

    Should an atom return itself as a neighbor?

  bothways: bool

    Return all neighbors. Default is to return only “half” of the neighbors.

  primitive: `PrimitiveNeighborList` or `NewPrimitiveNeighborList` class

    Define which implementation to use. Older and quadratically-scaling `PrimitiveNeighborList` or newer and linearly-scaling `NewPrimitiveNeighborList`.

  Example:

      nl = NeighborList([2.3, 1.7])
      nl.update(atoms)
      indices, offsets = nl.get_neighbors(0)

  get_connectivity_matrix(sparse=True)

    See `get_connectivity_matrix()`.

  get_neighbors(a)

    See `ase.neighborlist.PrimitiveNeighborList.get_neighbors()` or `ase.neighborlist.PrimitiveNeighborList.get_neighbors()`.

  property nneighbors

    Get number of neighbors.

  property npbcneighbors

    Get number of pbc neighbors.

  property nupdates

    Get number of updates.

  update(atoms)

    See `ase.neighborlist.PrimitiveNeighborList.update()` or `ase.neighborlist.PrimitiveNeighborList.update()`.

class ase.neighborlist.NewPrimitiveNeighborList(cutoffs, skin=0.3, sorted=False, self_interaction=True, bothways=False, use_scaled_positions=False)

  Neighbor list object. Wrapper around neighbor_list and first_neighbors.

  cutoffs: list of float

    List of cutoff radii - one for each atom. If the spheres (defined by their cutoff radii) of two atoms overlap, they will be counted as neighbors.

  skin: float

    If no atom has moved more than the skin-distance since the last call to the `update()` method, then the neighbor list can be reused. This will save some expensive rebuilds of the list, but extra neighbors outside the cutoff will be returned.

  sorted: bool

    Sort neighbor list.

  self_interaction: bool

    Should an atom return itself as a neighbor?

  bothways: bool

    Return all neighbors. Default is to return only “half” of the neighbors.

  Example:

      nl = NeighborList([2.3, 1.7])
      nl.update(atoms)
      indices, offsets = nl.get_neighbors(0)

  build(pbc, cell, positions, numbers=None)

    Build the list.

  get_neighbors(a)

    Return neighbors of atom number a.

    A list of indices and offsets to neighboring atoms is returned. The positions of the neighbor atoms can be calculated like this:

        >>>  indices, offsets = nl.get_neighbors(42)
        >>>  for i, offset in zip(indices, offsets):
        >>>      print(atoms.positions[i] + dot(offset, atoms.get_cell()))

    Notice that if get_neighbors(a) gives atom b as a neighbor, then get_neighbors(b) will not return a as a neighbor - unless bothways=True was used.

  update(pbc, cell, positions, numbers=None)

    Make sure the list is up to date.

class ase.neighborlist.PrimitiveNeighborList(cutoffs, skin=0.3, sorted=False, self_interaction=True, bothways=False, use_scaled_positions=False)

  Neighbor list that works without Atoms objects.

  This is less fancy, but can be used to avoid conversions between scaled and non-scaled coordinates which may affect cell offsets through rounding errors.

  build(pbc, cell, coordinates)

    Build the list.

    Coordinates are taken to be scaled or not according to self.use_scaled_positions.

  get_neighbors(a)

    Return neighbors of atom number a.

    A list of indices and offsets to neighboring atoms is returned. The positions of the neighbor atoms can be calculated like this:

        indices, offsets = nl.get_neighbors(42)
        for i, offset in zip(indices, offsets):
            print(atoms.positions[i] + offset @ atoms.get_cell())

    Notice that if get_neighbors(a) gives atom b as a neighbor, then get_neighbors(b) will not return a as a neighbor - unless bothways=True was used.

  update(pbc, cell, coordinates)

    Make sure the list is up to date.

ase.neighborlist.build_neighbor_list(atoms, cutoffs=None, **kwargs)

  Automatically build and update a NeighborList.

  Parameters:

  atoms`Atoms` object

    Atoms to build Neighborlist for.

  cutoffs: list of floats

    Radii for each atom. If not given it will be produced by calling `ase.neighborlist.natural_cutoffs()`

  kwargs: arbitrary number of options

    Will be passed to the constructor of `NeighborList`

  Returns:

  return: `NeighborList`

    A `NeighborList` instance (updated).

ase.neighborlist.first_neighbors(natoms, first_atom)

  Compute an index array pointing to the ranges within the neighbor list that contain the neighbors for a certain atom.

  Parameters:

  natomsint

    Total number of atom.

  first_atomarray_like

    Array containing the first atom ‘i’ of the neighbor tuple returned by the neighbor list.

  Returns:

  seedarray

    Array containing pointers to the start and end location of the neighbors of a certain atom. Neighbors of atom k have indices from s[k] to s[k+1]-1.

ase.neighborlist.get_connectivity_matrix(nl, sparse=True)

  Return connectivity matrix for a given NeighborList (dtype=numpy.int8).

  A matrix of shape (nAtoms, nAtoms) will be returned. Connected atoms i and j will have matrix[i,j] == 1, unconnected matrix[i,j] == 0. If bothways=True the matrix will be symmetric, otherwise not!

  If sparse is True, a scipy csr matrix is returned. If sparse is False, a numpy matrix is returned.

  Note that the old and new neighborlists might give different results for periodic systems if bothways=False.

  Example:

  Determine which molecule in a system atom 1 belongs to.

      >>> from ase import neighborlist
      >>> from ase.build import molecule
      >>> from scipy import sparse
      >>> mol = molecule('CH3CH2OH')
      >>> cutOff = neighborlist.natural_cutoffs(mol)
      >>> neighborList = neighborlist.NeighborList(cutOff, self_interaction=False, bothways=True)
      >>> neighborList.update(mol)
      >>> matrix = neighborList.get_connectivity_matrix()
      >>> #or: matrix = neighborlist.get_connectivity_matrix(neighborList.nl)
      >>> n_components, component_list = sparse.csgraph.connected_components(matrix)
      >>> idx = 1
      >>> molIdx = component_list[idx]
      >>> print("There are {} molecules in the system".format(n_components))
      >>> print("Atom {} is part of molecule {}".format(idx, molIdx))
      >>> molIdxs = [ i for i in range(len(component_list)) if component_list[i] == molIdx ]
      >>> print("The following atoms are part of molecule {}: {}".format(molIdx, molIdxs))

ase.neighborlist.get_distance_indices(distanceMatrix, distance)

  Get indices for each node that are distance or less away.

  Parameters:

  distanceMatrix: any one of scipy.sparse matrices (NxN)

    Matrix containing distance information of atoms. Get it e.g. with `get_distance_matrix()`

  distance: integer

    Number of steps to allow.

  Returns:

  return: list of length N

    A list of length N. return[i] has all indices that are connected to item i.

  The distance matrix only contains shortest paths, so when looking for distances longer than one, we need to add the lower values for cases where atoms are connected via a shorter path too.

ase.neighborlist.get_distance_matrix(graph, limit=3)

  Get Distance Matrix from a Graph.

  Parameters:

  graph: array, matrix or sparse matrix, 2 dimensions (N, N)

    Graph representation of the connectivity. See scipy doc for reference.

  limit: integer

    Maximum number of steps to analyze. For most molecular information, three should be enough.

  Returns:

  return: scipy.sparse.csr_matrix, shape (N, N)

    A scipy.sparse.csr_matrix. All elements that are not connected within limit steps are set to zero.

  This is a potential memory bottleneck, as csgraph.dijkstra produces a dense output matrix. Here we replace all np.inf values with 0 and transform back to csr_matrix. Why not dok_matrix like the connectivity-matrix? Because row-picking is most likely and this is super fast with csr.

ase.neighborlist.mic(dr, cell, pbc=True)

  Apply minimum image convention to an array of distance vectors.

  Parameters:

  drarray_like

    Array of distance vectors.

  cellarray_like

    Simulation cell.

  pbcarray_like, optional

    Periodic boundary conditions in x-, y- and z-direction. Default is to assume periodic boundaries in all directions.

  Returns:

  drarray

    Array of distance vectors, wrapped according to the minimum image convention.

ase.neighborlist.natural_cutoffs(atoms, mult=1, **kwargs)

  Generate a radial cutoff for every atom based on covalent radii.

  The covalent radii are a reasonable cutoff estimation for bonds in many applications such as neighborlists, so function generates an atoms length list of radii based on this idea.

  - atoms: An atoms object

  - mult: A multiplier for all cutoffs, useful for coarse grained adjustment

  - kwargs: Symbol of the atom and its corresponding cutoff, used to override the covalent radii

ase.neighborlist.neighbor_list(quantities, a, cutoff, self_interaction=False, max_nbins=1000000.0)

  Compute a neighbor list for an atomic configuration.

  Atoms outside periodic boundaries are mapped into the box. Atoms outside nonperiodic boundaries are included in the neighbor list but complexity of neighbor list search for those can become n^2.

  The neighbor list is sorted by first atom index ‘i’, but not by second atom index ‘j’.

  Parameters:

  quantities: str

    Quantities to compute by the neighbor list algorithm. Each character in this string defines a quantity. They are returned in a tuple of the same order. Possible quantities are:

    > - ‘i’ : first atom index

    > - ‘j’ : second atom index

    > - ‘d’ : absolute distance

    > - ‘D’ : distance vector

    > - ‘S’ : shift vector (number of cell boundaries crossed by the bond between atom i and j). With the shift vector S, the distances D between atoms can be computed from: D = a.positions[j]-a.positions[i]+S.dot(a.cell)

  a: `ase.Atoms`

    Atomic configuration.

  cutoff: float or dict

    Cutoff for neighbor search. It can be:

    > - A single float: This is a global cutoff for all elements.

    > - A dictionary: This specifies cutoff values for element pairs. Specification accepts element numbers of symbols. Example: {(1, 6): 1.1, (1, 1): 1.0, (‘C’, ‘C’): 1.85}

    > - A list/array with a per atom value: This specifies the radius of an atomic sphere for each atoms. If spheres overlap, atoms are within each others neighborhood. See `natural_cutoffs()` for an example on how to get such a list.

  self_interaction: bool

    Return the atom itself as its own neighbor if set to true. Default: False

  max_nbins: int

    Maximum number of bins used in neighbor search. This is used to limit the maximum amount of memory required by the neighbor list.

  Returns:

  i, j, …: array

    Tuple with arrays for each quantity specified above. Indices in \(i\) are returned in ascending order 0..len(a), but the order of (i,j) pairs is not guaranteed.

  Examples:

  Examples assume Atoms object a and numpy imported as np.

  1. Coordination counting:

         i = neighbor_list('i', a, 1.85)
         coord = np.bincount(i)

  2. Coordination counting with different cutoffs for each pair of species:

         i = neighbor_list('i', a,
                           {('H', 'H'): 1.1, ('C', 'H'): 1.3, ('C', 'C'): 1.85})
         coord = np.bincount(i)

  3. Pair distribution function:

         d = neighbor_list('d', a, 10.00)
         h, bin_edges = np.histogram(d, bins=100)
         pdf = h/(4*np.pi/3*(bin_edges[1:]**3 - bin_edges[:-1]**3)) * a.get_volume()/len(a)

  4. Pair potential:

         i, j, d, D = neighbor_list('ijdD', a, 5.0)
         energy = (-C/d**6).sum()
         pair_forces = (6*C/d**5  * (D/d).T).T
         forces_x = np.bincount(j, weights=pair_forces[:, 0], minlength=len(a)) -                    np.bincount(i, weights=pair_forces[:, 0], minlength=len(a))
         forces_y = np.bincount(j, weights=pair_forces[:, 1], minlength=len(a)) -                    np.bincount(i, weights=pair_forces[:, 1], minlength=len(a))
         forces_z = np.bincount(j, weights=pair_forces[:, 2], minlength=len(a)) -                    np.bincount(i, weights=pair_forces[:, 2], minlength=len(a))

  5. Dynamical matrix for a pair potential stored in a block sparse format:

         from scipy.sparse import bsr_matrix
         i, j, dr, abs_dr = neighbor_list('ijDd', atoms)
         energy = (dr.T / abs_dr).T
         dynmat = -(dde * (energy.reshape(-1, 3, 1) * energy.reshape(-1, 1, 3)).T).T                  -(de / abs_dr * (np.eye(3, dtype=energy.dtype) -                    (energy.reshape(-1, 3, 1) * energy.reshape(-1, 1, 3))).T).T
         dynmat_bsr = bsr_matrix((dynmat, j, first_i), shape=(3*len(a), 3*len(a)))

         dynmat_diag = np.empty((len(a), 3, 3))
         for x in range(3):
             for y in range(3):
                 dynmat_diag[:, x, y] = -np.bincount(i, weights=dynmat[:, x, y])

         dynmat_bsr += bsr_matrix((dynmat_diag, np.arange(len(a)),
                                   np.arange(len(a) + 1)),
                                  shape=(3 * len(a), 3 * len(a)))

ase.neighborlist.primitive_neighbor_list(quantities, pbc, cell, positions, cutoff, numbers=None, self_interaction=False, use_scaled_positions=False, max_nbins=1000000.0)

  Compute a neighbor list for an atomic configuration.

  Atoms outside periodic boundaries are mapped into the box. Atoms outside nonperiodic boundaries are included in the neighbor list but complexity of neighbor list search for those can become n^2.

  The neighbor list is sorted by first atom index ‘i’, but not by second atom index ‘j’.

  Parameters:

  quantities: str

    Quantities to compute by the neighbor list algorithm. Each character in this string defines a quantity. They are returned in a tuple of the same order. Possible quantities are

    > - ‘i’ : first atom index

    > - ‘j’ : second atom index

    > - ‘d’ : absolute distance

    > - ‘D’ : distance vector

    > - ‘S’ : shift vector (number of cell boundaries crossed by the bond between atom i and j). With the shift vector S, the distances D between atoms can be computed from: D = positions[j]-positions[i]+S.dot(cell)

  pbc: array_like

    3-tuple indicating giving periodic boundaries in the three Cartesian directions.

  cell: 3x3 matrix

    Unit cell vectors.

  positions: list of xyz-positions

    Atomic positions. Anything that can be converted to an ndarray of shape (n, 3) will do: [(x1,y1,z1), (x2,y2,z2), …]. If use_scaled_positions is set to true, this must be scaled positions.

  cutoff: float or dict

    Cutoff for neighbor search. It can be:

    > - A single float: This is a global cutoff for all elements.

    > - A dictionary: This specifies cutoff values for element pairs. Specification accepts element numbers of symbols. Example: {(1, 6): 1.1, (1, 1): 1.0, (‘C’, ‘C’): 1.85}

    > - A list/array with a per atom value: This specifies the radius of an atomic sphere for each atoms. If spheres overlap, atoms are within each others neighborhood. See `natural_cutoffs()` for an example on how to get such a list.

  self_interaction: bool

    Return the atom itself as its own neighbor if set to true. Default: False

  use_scaled_positions: bool

    If set to true, positions are expected to be scaled positions.

  max_nbins: int

    Maximum number of bins used in neighbor search. This is used to limit the maximum amount of memory required by the neighbor list.

  Returns:

  i, j, …array

    Tuple with arrays for each quantity specified above. Indices in \(i\) are returned in ascending order 0..len(a)-1, but the order of (i,j) pairs is not guaranteed.