Neighbour lists

class matscipy.neighbours.Neighbourhood(atom_types=None)

Bases: ABC

Abstract class defining a neighbourhood of atoms (pairs, triplets).

__init__(atom_types=None)

Initialize with atoms and optional atom types.

abstract get_pairs(atoms: Atoms, quantities: str, cutoff=None)

Return requested data on pairs.

abstract get_triplets(atoms: Atoms, quantities: str, neighbours=None, cutoff=None, full_connectivity=False)

Return requested data on triplets.

static mask(pair_distances, cutoff)
static make_result(quantities, connectivity, D, d, S, accepted_quantities) → List

Construct result list.

static compute_distances(atoms: Atoms, connectivity: ndarray, indices: List[int]) → Tuple[ndarray, ndarray]

Return distances and vectors for connectivity.

connected_triplets(atoms: Atoms, pair_list, triplet_list, nb_pairs)
triplet_to_numbers(atoms: Atoms, i, j, k)
find_triplet_types(atoms: Atoms, i, j, k)

Return triplet types from atom ids.

static lexsort(connectivity: ndarray)
abstract double_neighbourhood()

Return neighbourhood with double cutoff/connectivity.

abstract reverse_pair_indices(i_p: ndarray, j_p: ndarray, r_p: ndarray)

Return indices of reverse pairs.

class matscipy.neighbours.CutoffNeighbourhood(atom_types=None, pair_types=None, triplet_types=None, cutoff: float | dict | None = None)

Bases: Neighbourhood

Class defining neighbourhood based on proximity.

__init__(atom_types=None, pair_types=None, triplet_types=None, cutoff: float | dict | None = None)

Initialize with atoms, atom types, pair types and cutoff.

Parameters:

  • atom_types (ArrayLike) – atom types array

  • pair_types (function of 2 atom type arrays) – maps 2 atom types array to an array of pair types

  • 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.

get_pairs(atoms: Atoms, quantities: str, cutoff=None)

Return pairs and quantities from conventional neighbour list.

get_triplets(atoms: Atoms, quantities: str, neighbours=None, cutoff=None)

Return triplets and quantities from conventional neighbour list.

double_neighbourhood()

Return neighbourhood with double cutoff/connectivity.

reverse_pair_indices(i_p: ndarray, j_p: ndarray, r_p: ndarray)

Return indices of reverse pairs.

static compute_distances(atoms: Atoms, connectivity: ndarray, indices: List[int]) → Tuple[ndarray, ndarray]

Return distances and vectors for connectivity.

connected_triplets(atoms: Atoms, pair_list, triplet_list, nb_pairs)
find_triplet_types(atoms: Atoms, i, j, k)

Return triplet types from atom ids.

static lexsort(connectivity: ndarray)
static make_result(quantities, connectivity, D, d, S, accepted_quantities) → List

Construct result list.

static mask(pair_distances, cutoff)
triplet_to_numbers(atoms: Atoms, i, j, k)
class matscipy.neighbours.MolecularNeighbourhood(molecules: Molecules, atom_types=None, double_cutoff=False)

Bases: Neighbourhood

Class defining neighbourhood based on molecular connectivity.

__init__(molecules: Molecules, atom_types=None, double_cutoff=False)

Initialze with atoms and molecules.

double_neighbourhood()

Return neighbourhood with double cutoff/connectivity.

property molecules

Molecules instance that defines neighbourhood.

property pair_type

Map atom types to pair types.

property triplet_type

Map atom types to triplet types.

static double_connectivity(connectivity: ndarray) → ndarray

Sort and stack connectivity + reverse connectivity.

complete_connectivity(typeoffset: int = 0)

Add angles to pair connectivity.

get_pairs(atoms: Atoms, quantities: str, cutoff=None)

Return pairs and quantities from connectivities.

get_triplets(atoms: Atoms, quantities: str, neighbours=None, cutoff=None)

Return triplets and quantities from connectivities.

find_triplet_types(atoms: Atoms, i, j, k)

Return triplet types from atom ids.

reverse_pair_indices(i_p: ndarray, j_p: ndarray, r_p: ndarray)

Return indices of reverse pairs.

static compute_distances(atoms: Atoms, connectivity: ndarray, indices: List[int]) → Tuple[ndarray, ndarray]

Return distances and vectors for connectivity.

connected_triplets(atoms: Atoms, pair_list, triplet_list, nb_pairs)
static lexsort(connectivity: ndarray)
static make_result(quantities, connectivity, D, d, S, accepted_quantities) → List

Construct result list.

static mask(pair_distances, cutoff)
triplet_to_numbers(atoms: Atoms, i, j, k)
matscipy.neighbours.mic(dr, cell, pbc=None)

Apply minimum image convention to an array of distance vectors.

Parameters:

  • dr (array_like) – Array of distance vectors.

  • cell (array_like) – Simulation cell.

  • pbc (array_like, optional) – Periodic boundary conditions in x-, y- and z-direction. Default is to assume periodic boundaries in all directions.

Returns:

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

Return type:

array

matscipy.neighbours.neighbour_list(quantities, atoms=None, cutoff=None, positions=None, cell=None, pbc=None, numbers=None, cell_origin=None)

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 the 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’.

The neighbour list accepts either an ASE Atoms object or positions and cell vectors individually.

quantitiesstr

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)

atomsase.Atoms

Atomic configuration. (Default: None)

cutofffloat 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.

positionsarray_like

Atomic positions. (Default: None)

cellarray_like

Cell vectors as a 3x3 matrix. (Default: Shrink wrapped cell)

pbcarray_like

3-vector containing periodic boundary conditions in all three directions. (Default: Nonperiodic box)

numbersarray_like

Array containing the atomic numbers.

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 assume Atoms object a and numpy imported as np. 1. Coordination counting:

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

  1. 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)

  2. 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)

  3. 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))

  4. 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)))

i_n, j_n, dr_nc, abs_dr_n = neighbour_list(‘ijDd’, atoms, dict)

e_nc = (dr_nc.T/abs_dr_n).T

D_ncc = -(dde_n * (e_nc.reshape(-1,3,1) * e_nc.reshape(-1,1,3)).T).T D_ncc += -(de_n/abs_dr_n * (np.eye(3, dtype=e_nc.dtype) - (e_nc.reshape(-1,3,1) * e_nc.reshape(-1,1,3))).T).T

D = bsr_matrix((D_ncc, j_n, first_i), shape=(3*nat,3*nat))

Ddiag_icc = np.empty((nat,3,3)) for x in range(3):

for y in range(3):

Ddiag_icc[:,x,y] = -np.bincount(i_n, weights = D_ncc[:,x,y])

D += bsr_matrix((Ddiag_icc,np.arange(nat),np.arange(nat+1)), shape=(3*nat,3*nat))

return D

matscipy.neighbours.triplet_list(first_neighbours, abs_dr_p=None, cutoff=None, i_p=None, j_p=None)

Compute a triplet list for an atomic configuration. The triple list is a mask that can be applied to the corresponding neighbour list to mask triplet properties. The triplet list accepts an first_neighbour array (generated by first_neighbours) as input.

Parameters:

first_neighbours (array) – adresses of the first time an atom occours in the neighour list

Returns:

  • ij_t, ik_t (array) – lists of adresses that form triples in the pair lists

  • jk_t (array (if and only if i_p, j_p, first_i != None)) – list of pairs jk that connect each triplet ij, ik between atom j and k

Example

i_n, j_n, abs_dr_p = neighbour_list(‘ijd’, atoms=atoms, cutoff=cutoff)

first_i = np.array([0, 2, 6, 10], dtype=’int32’) a = triplet_list(first_i, [2.2]*9+[3.0], 2.6)

# one may obtain first_ij by using find_ij = first_neighbours(len(i_p), ij_t) # or (slower but less parameters and more general, # i.e for every ordered list) first_ij = get_jump_indicies(ij_t)

matscipy.neighbours.find_indices_of_reversed_pairs(i_n, j_n, abs_dr_n)

Find neighbor list indices where reversed pairs are stored

Given list of identifiers of neighbor atoms i_n and j_n, determines the list of indices reverse into the neighbor list where each pair is reversed, i.e. i_n[reverse[n]]=j_n[n] and j_n[reverse[n]]=i_n[n] for each index n in the neighbor list

In the case of small periodic systems, one needs to be careful, because the same pair may appear more than one time, with different pair distances. Therefore, the pair distance must be taken into account.

We assume that there is in fact one reversed pair for every pair. However, we do not check this assumption in order to avoid overhead.

Parameters:

  • i_n (array_like) – array of atom identifiers

  • j_n (array_like) – array of atom identifiers

  • abs_dr_n (array_like) – pair distances

Returns:

reverse – array of indices into i_n and j_n

Return type:

numpy.ndarray

matscipy.neighbours.find_common_neighbours(i_n, j_n, nat)

Find common neighbors of pairs of atoms

For each pair (i1, j1) in the neighbor list, find all other pairs (i2, j1) which share the same j1. This includes (i1,j1) itself. In this way, create a list with n blocks of rows, where n is the length of the neighbor list. All rows in a block have the same j1. Each row corresponds to one triplet (i1, j1 ,i2). The number of rows in the block is equal to the total number of neighbors of j1.

Parameters:

  • i_n (array_like) – array of atom identifiers

  • j_n (array_like) – array of atom identifiers

  • nat (int) – number of atoms

Returns:

  • cnl_i1_i2 (array) – atom numbers i1 and i2

  • cnl_j1 (array) – shared neighbor of i1 and i2

  • nl_index_i1_j1 (array) – index in the neighbor list of pair i1, j1

  • nl_index_i2_j1 (array) – index in the neighbor list of pair i2, j1

Examples

Accumulate random numbers for pairs with common neighbors:

>>> import numpy as np
>>> import matscipy
>>> from ase.lattice.cubic import FaceCenteredCubic
>>> from matscipy.neighbours import neighbour_list, find_common_neighbours
>>> cutoff = 6.5
>>> atoms = FaceCenteredCubic('Cu', size=[4, 4, 4])
>>> nat = len(atoms.numbers)
>>> print(nat)
256
>>> i_n, j_n, dr_nc, abs_dr_n = neighbour_list('ijDd', atoms, cutoff)
>>> print(i_n.shape)
(22016,)
>>> cnl_i1_i2, cnl_j1, nl_index_i1_j1, nl_index_i2_j1 = find_common_neighbours(i_n, j_n, nat)
>>> print(cnl_i1_i2.shape)
(1893376, 2)
>>> unique_pairs_i1_i2, bincount_bins = np.unique(cnl_i1_i2, axis=0, return_inverse=True)
>>> print(unique_pairs_i1_i2.shape)
(65536, 2)
>>> tmp = np.random.rand(cnl_i1_i2.shape[0])
>>> my_sum = np.bincount(bincount_bins, weights=tmp, minlength=unique_pairs_i1_i2.shape[0])
>>> print(my_sum.shape)
(65536,)