Neighbour list

Basic usage

matscipy neighbour lists are stored in a format similar to the coordinate (COO) format of sparse matrices. The basic neighbor list consists of two array that each contain the indices of the atoms that constitute the pair.

import numpy as np

from ase.build import molecule
from matscipy.neighbours import neighbour_list

# Single water in a box with vacuum
a = molecule('H2O')
a.center(vacuum=5)

# Create neighbor list
i, j = neighbour_list('ij', a, cutoff=1.2)

# Return list of neighbor pairs
a.get_chemical_symbols(), i, j

(['O', 'H', 'H'],
 array([0, 0, 1, 2], dtype=int32),
 array([2, 1, 0, 0], dtype=int32))

The water molecule has four pairs at a cutoff of 1.2, which are the O-H bonds. Each of the bonds occurs twice in the neighbor list.

This list format allows simple analysis. For example, coordination numbers can be computed by counting the number of entries in the index arrays.

# Bincount counts the number of times a specific entry shows up in the array
np.bincount(i)

array([2, 1, 1])

The oxygen atom has a coordination of 2 (both hydrogens) while each of the hydrogens has a coordination of 1 (since only the oxygen is the neighbor).

Per-atom properties

The neighbour list can also compute per atom properties, in particular distances and distance vectors. The first argument to the neighbour_list function is a string that identifies the members of the return tuple. If we want distances between atoms, we additionally specific a ‘d’ in this string. The return tuple then has three members.

i, j, d = neighbour_list('ijd', a, cutoff=1.2)

d

array([0.96856502, 0.96856502, 0.96856502, 0.96856502])

This the O-H bond length. If we increase the cutoff to 2 Å, we also capture the H-H distance.

neighbour_list('d', a, cutoff=2.0)

array([0.96856502, 0.96856502, 0.96856502, 1.526478  , 0.96856502,
       1.526478  ])

Interatomic potential

As an advances usage of the neighbor list, consider the implementation of a pair potential on top of the neighbour list data structure. (This is actually how it is done in the calculators that ship with matscipy.) The following code example implements an attractive \(\propto r^{-6}\) potential, i.e. a London dispersion force with a prefactor {math}:C.

from matscipy.numpy_tricks import mabincount

C = 1.0  # just some number

i, j, d, D = neighbour_list('ijdD', a, 5.0)
energy = (-C/d**6).sum()
pair_forces = (6*C/d**5  * D.T/d).T
forces = mabincount(j, pair_forces, len(a)) - mabincount(i, pair_forces, len(a))

forces

array([[  0.        ,   0.        ,  17.33444449],
       [  0.        ,  12.54138009,  -8.66722225],
       [  0.        , -12.54138009,  -8.66722225]])

The code computes the energy as a sum over all pair contributions. The variable pair_forces contains the force vectors between pairs of atoms, which are then summed onto the respective components of the force. The utility function mabincount works like np.bincount but can handle multidimensional arrays.