cna

Common neighbour analysis (CNA) of local crystal structure.

Each bond from an atom to one of its neighbours is labelled by three numbers: how many neighbours the pair has in common (n_cn), how many bonds those common neighbours share (n_b), and how many bonds the largest connected cluster of those holds (n_lcb). An atom's multiset of labels identifies its local structure:

FCC  12 x (4, 2, 1)
HCP   6 x (4, 2, 1) + 6 x (4, 2, 2)
ICO  12 x (5, 5, 5)
BCC   8 x (6, 6, 6) + 6 x (4, 4, 4)

Everything about one atom happens inside its own neighbour shell, so the analysis is built from layers, each usable on its own:

• chmpy.calc.nearest.nearest_neighbors finds the shell.
• neighbor_bonds records which neighbours are bonded to which as bitmasks, one integer per neighbour, bit b of entry a set when neighbours a and b are bonded. The common neighbours of the atom and neighbour a are then exactly the set bits of entry a.
• cna_signatures turns the bitmasks into (n_cn, n_b, n_lcb), and classify matches them against PATTERNS.

Three ways of drawing the bonds are provided:

• fixed_cna uses one cutoff throughout (Honeycutt and Andersen, J. Phys. Chem. 91, 4950 (1987)).
• adaptive_cna sets each atom's cutoff from its own neighbour distances (Stukowski, Modelling Simul. Mater. Sci. Eng. 20, 045021 (2012)).
• interval_cna needs no cutoff, and tolerates far more thermal disorder (Larsen, arXiv:2003.08879). Bonds of equal length do not bound an interval (INTERVAL_TOLERANCE), so that the answer for an atom does not depend on the rounding error between them, or on the orientation of the structure.

Atoms are processed in chunks, on workers threads; numpy releases the GIL for the work inside each chunk.

StructureType

Bases: IntEnum

Local structures CNA identifies, with their conventional numbering.

Source code in chmpy/descriptors/cna.py

class StructureType(IntEnum):
    """Local structures CNA identifies, with their conventional numbering."""

    OTHER = 0
    FCC = 1
    HCP = 2
    BCC = 3
    ICO = 4

adaptive_cna(positions, cell=None, pbc=None, neighbors=None, workers=-1)

Adaptive common neighbour analysis.

FCC, HCP and icosahedral atoms are tested on their 12 nearest neighbours with a cutoff of (1 + sqrt 2) / 2 times their mean distance; atoms that are none of those are tested for BCC on their 14 nearest, with the cutoff scaled from the lattice parameter those imply.

Parameters:

Name	Type	Description	Default
positions		(N, 3) Cartesian positions, Angstroms	required
cell		(3, 3) lattice vectors as rows, or None for an isolated system	None
pbc		(3,) periodic directions, as for nearest_neighbors	None
neighbors		NearestNeighbors with at least 14 per atom, to reuse a search already done; found from the other arguments if omitted	None
workers		threads to use, -1 for one per CPU	-1

Returns:

Type	Description
	(N,) int8 StructureType of each atom

Source code in chmpy/descriptors/cna.py

def adaptive_cna(positions, cell=None, pbc=None, neighbors=None, workers=-1):
    """Adaptive common neighbour analysis.

    FCC, HCP and icosahedral atoms are tested on their 12 nearest neighbours
    with a cutoff of `(1 + sqrt 2) / 2` times their mean distance; atoms that
    are none of those are tested for BCC on their 14 nearest, with the cutoff
    scaled from the lattice parameter those imply.

    Args:
        positions: (N, 3) Cartesian positions, Angstroms
        cell: (3, 3) lattice vectors as rows, or None for an isolated system
        pbc: (3,) periodic directions, as for `nearest_neighbors`
        neighbors: `NearestNeighbors` with at least 14 per atom, to reuse a
            search already done; found from the other arguments if omitted
        workers: threads to use, -1 for one per CPU

    Returns:
        (N,) int8 `StructureType` of each atom
    """
    neighbors = _neighbors(neighbors, 14, positions, cell, pbc, workers, "adaptive")
    vectors, distances = neighbors.vectors, neighbors.distances
    types = np.zeros(len(neighbors), dtype=np.int8)

    atoms = np.flatnonzero(np.isfinite(distances[:, 11]))
    cutoff = distances[atoms, :12].mean(axis=1) * FCC_CUTOFF_FACTOR
    bonds = neighbor_bonds(vectors[atoms, :12], cutoff, workers)
    types[atoms] = classify(bonds, PATTERNS[12], workers)

    atoms = np.flatnonzero(
        (types == StructureType.OTHER) & np.isfinite(distances[:, 13])
    )
    first, second = distances[atoms, :8], distances[atoms, 8:14]
    lattice = (first.sum(axis=1) / np.sqrt(0.75) + second.sum(axis=1)) / 14
    bonds = neighbor_bonds(vectors[atoms, :14], lattice * BCC_CUTOFF_FACTOR, workers)
    types[atoms] = classify(bonds, PATTERNS[14], workers)
    return types

classify(bonds, patterns, workers=-1)

Match each atom's signatures against a set of structures.

Only atoms whose common-neighbour counts fit some structure have their signatures computed, and only the bonds whose first two numbers fit have their chains followed, so a disordered atom costs little.

Parameters:

Name	Type	Description	Default
bonds		(N, k) bitmasks from neighbor_bonds	required
patterns		{structure: {(n_cn, n_b, n_lcb): count}}, the counts of each structure summing to k, as in PATTERNS	required
workers		threads to use, -1 for one per CPU	-1

Returns:

Type	Description
ndarray	(N,) int8 structure of each atom, StructureType.OTHER where none
ndarray	matches

Source code in chmpy/descriptors/cna.py

def classify(bonds, patterns, workers=-1) -> np.ndarray:
    """Match each atom's signatures against a set of structures.

    Only atoms whose common-neighbour counts fit some structure have their
    signatures computed, and only the bonds whose first two numbers fit have
    their chains followed, so a disordered atom costs little.

    Args:
        bonds: (N, k) bitmasks from `neighbor_bonds`
        patterns: `{structure: {(n_cn, n_b, n_lcb): count}}`, the counts of
            each structure summing to k, as in `PATTERNS`
        workers: threads to use, -1 for one per CPU

    Returns:
        (N,) int8 structure of each atom, `StructureType.OTHER` where none
        matches
    """
    bonds = np.asarray(bonds, dtype=np.uint32)
    out = np.zeros(len(bonds), dtype=np.int8)

    def work(start, stop):
        out[start:stop] = _classify_chunk(bonds[start:stop], patterns)

    _for_chunks(len(bonds), work, workers)
    return out

cna_signatures(bonds, workers=-1)

The CNA signature of every bond from every atom.

Parameters:

Name	Type	Description	Default
bonds		(N, k) bitmasks from neighbor_bonds	required
workers		threads to use, -1 for one per CPU	-1

Returns:

Type	Description
ndarray	(N, k, 3) int (n_cn, n_b, n_lcb) for the bond from each atom to each
ndarray	of its neighbours

Source code in chmpy/descriptors/cna.py

def cna_signatures(bonds, workers=-1) -> np.ndarray:
    """The CNA signature of every bond from every atom.

    Args:
        bonds: (N, k) bitmasks from `neighbor_bonds`
        workers: threads to use, -1 for one per CPU

    Returns:
        (N, k, 3) int `(n_cn, n_b, n_lcb)` for the bond from each atom to each
        of its neighbours
    """
    bonds = np.asarray(bonds, dtype=np.uint32)
    n, k = bonds.shape
    out = np.empty((n, k, 3), dtype=np.int64)

    def work(start, stop):
        chunk = bonds[start:stop]
        common, shared, chain = _signature(chunk)
        _resolve_chains(chunk, chain)
        out[start:stop] = np.stack([common, shared, chain], axis=-1)

    _for_chunks(n, work, workers)
    return out

fixed_cna(positions, cutoff, cell=None, pbc=None, neighbors=None, workers=-1)

Conventional common neighbour analysis with one cutoff.

Atoms with exactly 12 neighbours inside the cutoff are tested for FCC, HCP and icosahedral order, those with exactly 14 for BCC, and the rest are StructureType.OTHER.

Parameters:

Name	Type	Description	Default
positions		(N, 3) Cartesian positions, Angstroms	required
cutoff		bond length cutoff, Angstroms	required
cell		(3, 3) lattice vectors as rows, or None for an isolated system	None
pbc		(3,) periodic directions, as for nearest_neighbors	None
neighbors		NearestNeighbors with at least 15 per atom, to reuse a search already done; found from the other arguments if omitted	None
workers		threads to use, -1 for one per CPU	-1

Returns:

Type	Description
	(N,) int8 StructureType of each atom

Source code in chmpy/descriptors/cna.py

def fixed_cna(positions, cutoff, cell=None, pbc=None, neighbors=None, workers=-1):
    """Conventional common neighbour analysis with one cutoff.

    Atoms with exactly 12 neighbours inside the cutoff are tested for FCC,
    HCP and icosahedral order, those with exactly 14 for BCC, and the rest are
    `StructureType.OTHER`.

    Args:
        positions: (N, 3) Cartesian positions, Angstroms
        cutoff: bond length cutoff, Angstroms
        cell: (3, 3) lattice vectors as rows, or None for an isolated system
        pbc: (3,) periodic directions, as for `nearest_neighbors`
        neighbors: `NearestNeighbors` with at least 15 per atom, to reuse a
            search already done; found from the other arguments if omitted
        workers: threads to use, -1 for one per CPU

    Returns:
        (N,) int8 `StructureType` of each atom
    """
    neighbors = _neighbors(neighbors, 15, positions, cell, pbc, workers, "fixed")
    coordination = (neighbors.distances <= cutoff).sum(axis=1)
    types = np.zeros(len(neighbors), dtype=np.int8)
    for shell, patterns in PATTERNS.items():
        atoms = np.flatnonzero(coordination == shell)
        bonds = neighbor_bonds(neighbors.vectors[atoms, :shell], cutoff, workers)
        types[atoms] = classify(bonds, patterns, workers)
    return types

interval_cna(positions, cell=None, pbc=None, neighbors=None, workers=-1)

Interval common neighbour analysis, which needs no cutoff.

Rather than fix a cutoff, the bonds among an atom's neighbours are added shortest first. Each structure is recognised at the point where the bonds drawn so far give its signatures, and the gap to the next bond is the interval of cutoffs over which that holds. The structure with the widest interval wins.

The 12 nearest neighbours are tested for FCC, HCP and icosahedral order and the 14 nearest for BCC. Only bonds shorter than two thirds of the way from the first neighbour shell to the second are considered, and the interval after the last of those ends at a longer bond chosen by _interval_ends.

Parameters:

Name	Type	Description	Default
positions		(N, 3) Cartesian positions, Angstroms	required
cell		(3, 3) lattice vectors as rows, or None for an isolated system	None
pbc		(3,) periodic directions, as for nearest_neighbors	None
neighbors		NearestNeighbors with at least 14 per atom, to reuse a search already done; found from the other arguments if omitted	None
workers		threads to use, -1 for one per CPU	-1

Returns:

Type	Description
	(N,) int8 StructureType of each atom

Source code in chmpy/descriptors/cna.py

def interval_cna(positions, cell=None, pbc=None, neighbors=None, workers=-1):
    """Interval common neighbour analysis, which needs no cutoff.

    Rather than fix a cutoff, the bonds among an atom's neighbours are added
    shortest first. Each structure is recognised at the point where the bonds
    drawn so far give its signatures, and the gap to the next bond is the
    interval of cutoffs over which that holds. The structure with the widest
    interval wins.

    The 12 nearest neighbours are tested for FCC, HCP and icosahedral order
    and the 14 nearest for BCC. Only bonds shorter than two thirds of the way
    from the first neighbour shell to the second are considered, and the
    interval after the last of those ends at a longer bond chosen by
    `_interval_ends`.

    Args:
        positions: (N, 3) Cartesian positions, Angstroms
        cell: (3, 3) lattice vectors as rows, or None for an isolated system
        pbc: (3,) periodic directions, as for `nearest_neighbors`
        neighbors: `NearestNeighbors` with at least 14 per atom, to reuse a
            search already done; found from the other arguments if omitted
        workers: threads to use, -1 for one per CPU

    Returns:
        (N,) int8 `StructureType` of each atom
    """
    neighbors = _neighbors(neighbors, 14, positions, cell, pbc, workers, "interval")
    distances = neighbors.distances
    types = np.zeros(len(neighbors), dtype=np.int8)
    full = np.isfinite(distances[:, 13])
    for shell, atoms in (
        (14, np.flatnonzero(full)),
        (12, np.flatnonzero(~full & np.isfinite(distances[:, 11]))),
    ):
        flat, n, _ = _flatten(neighbors.vectors[atoms, :shell])
        lengths = distances[atoms, :shell]
        out = np.zeros(n, dtype=np.int8)

        def work(start, stop, flat=flat, lengths=lengths, out=out):
            out[start:stop] = _interval_chunk(flat[start:stop], lengths[start:stop])

        _for_chunks(n, work, workers)
        types[atoms] = out
    return types

neighbor_bonds(vectors, cutoff, workers=-1)

Which neighbours of each atom are bonded to each other.

Parameters:

Name	Type	Description	Default
vectors		(N, k) neighbour vectors from each atom, as in NearestNeighbors.vectors; NaN entries bond to nothing	required
cutoff		bond length cutoff, a scalar or one per atom (N,)	required
workers		threads to use, -1 for one per CPU	-1

Returns:

Type	Description
ndarray	(N, k) uint32 bitmasks: bit b of entry [i, a] is set when
ndarray	neighbours a and b of atom i are no further apart than the cutoff

Source code in chmpy/descriptors/cna.py

def neighbor_bonds(vectors, cutoff, workers=-1) -> np.ndarray:
    """Which neighbours of each atom are bonded to each other.

    Args:
        vectors: (N, k) neighbour vectors from each atom, as in
            `NearestNeighbors.vectors`; NaN entries bond to nothing
        cutoff: bond length cutoff, a scalar or one per atom (N,)
        workers: threads to use, -1 for one per CPU

    Returns:
        (N, k) uint32 bitmasks: bit `b` of entry `[i, a]` is set when
        neighbours `a` and `b` of atom `i` are no further apart than the cutoff
    """
    flat, n, k = _flatten(vectors)
    cutoff2 = np.broadcast_to(np.asarray(cutoff, dtype=np.float64) ** 2, (n,))
    out = np.empty((n, k), dtype=np.uint32)

    def work(start, stop):
        distance2 = _pair_distances2(flat[start:stop], k)
        out[start:stop] = _pack(distance2 <= cutoff2[start:stop, None], k)

    _for_chunks(n, work, workers)
    return out