Coordination

This module contains functions to compute the coordination number of atoms in a system, including generalized coordination numbers and their strained editions.

snow.descriptors.coordination

coordination_number(coords, cut_off, neigh_list=False, pbc=False, box=None)

Computes the coordination number (number of nearest neighbours within a cutoff) for each atom in the system.

Optionally, the neighbour list (which is the heavy part of the calculation) can be returned with neigh_list=True.

Parameters:

Name	Type	Description	Default
coords	ndarray	Array of the coordinates of the atoms forming the system.	required
cut_off	float	The cutoff distance for determining nearest neighbors.	required
neigh_list	bool	Option to return the neighbour list as well as the coordination number of each atom	False
pbc	bool	Whether to apply periodic boundary conditions.	False
box	ndarray	Simulation box size (either [Lx, Ly, Lz] or [[xmin, xmax], [ymin, ymax], [zmin, zmax]] or 3 cell vectors (shape (3,3) - slower)).	None

Returns:

Name	Type	Description
neighs	list	neighbour list, the list of indeces of the neighbours of each atom. Only returned if neigh_list is True
coord_num	ndarray	the coordination numbers of each atom

Source code in snow/descriptors/coordination.py

def coordination_number(coords, cut_off, neigh_list=False, pbc=False, box=None):
    """
    Computes the coordination number (number of nearest neighbours within a cutoff) for each atom in the system.

    Optionally, the neighbour list (which is the heavy part of the calculation) can be returned with neigh_list=True.


    Parameters
    ----------
    coords : np.ndarray
        Array of the coordinates of the atoms forming the system.
    cut_off : float
        The cutoff distance for determining nearest neighbors.
    neigh_list : bool, default False
        Option to return the neighbour list as well as the coordination number of each atom
    pbc : bool, default False
        Whether to apply periodic boundary conditions.
    box : np.ndarray, optional
        Simulation box size (either [Lx, Ly, Lz] or [[xmin, xmax], [ymin, ymax], [zmin, zmax]] or 3 cell vectors (shape (3,3) - slower)).

    Returns
    -------
    neighs : list
        neighbour list, the list of indeces of the neighbours of each atom. Only returned if `neigh_list` is True
    coord_num : ndarray
        the coordination numbers of each atom
    """
    neigh = nearest_neighbours(
        coords=coords, cut_off=cut_off, pbc=pbc, box=box
    )
    n_atoms = np.shape(coords)[0]
    coord_numb = np.zeros(n_atoms)
    for i in range(n_atoms):
        coord_numb[i] = len(neigh[i])
    if neigh_list:
        return neigh, coord_numb
    else:
        return coord_numb

agcn_calculator(coords, cut_off, cn_max=12.0, strained=False, pbc=False, box=None, **kwargs)

Calculates the atop Generalized Coordination Number (GCN) for a site. The GCN is defined as the sum of the coordination numbers of the neighbors of each atom divided by the maximum typical coordination number in the specific system (cn_max).

Parameters:

Name	Type	Description	Default
coords	ndarray	Array of the coordinates of the atoms forming the system.	required
cut_off	float	The cutoff distance for determining nearest neighbors.	required
cn_max	float	Maximum coordination number in the specific system (default is 12.0, ok for fcc materials).	12.0
strained	bool	if True, computes the strained aGCN.	False
pbc	bool	Whether to apply periodic boundary conditions. Defaults to False.	False
box	ndarray	Simulation box size (either [Lx, Ly, Lz] or [[xmin, xmax], [ymin, ymax], [zmin, zmax]] or 3 cell vectors (shape (3,3) - slower)).	None
kwargs		thr_cn: int, optional a threshold coordination number value. If provided, only atoms with coordination < thr_cn are considered for the GCN calculation (e.g. useful if you only want to consider surface atoms in your calculation) dbulk: float, optional Bulk distance for strained aGCN (default is None), has to be provided if strained is True	{}

Returns:

Name	Type	Description
sites	ndarray	sites (atomic coordinates. If thr_cn is enabled, only coordinates of atoms selected by the coordination condition)
gcns	ndarray	Values of the atop GCN.

Source code in snow/descriptors/coordination.py

def agcn_calculator(coords, cut_off, cn_max = 12.0, strained: bool = False, pbc: bool = False, box = None, **kwargs):
    """
    Calculates the atop Generalized Coordination Number (GCN) for a site. The GCN is defined as the sum of the coordination numbers of the neighbors
    of each atom divided by the maximum typical coordination number in the specific system (cn_max).

    Parameters
    ----------
    coords : ndarray
        Array of the coordinates of the atoms forming the system.
    cut_off : float
        The cutoff distance for determining nearest neighbors.
    cn_max : float, default 12.0
        Maximum coordination number in the specific system (default is 12.0, ok for fcc materials).
    strained : bool, default False
        if True, computes the strained aGCN.
    pbc : bool, default False
        Whether to apply periodic boundary conditions. Defaults to False.
    box : np.ndarray, optional
        Simulation box size (either [Lx, Ly, Lz] or [[xmin, xmax], [ymin, ymax], [zmin, zmax]] or 3 cell vectors (shape (3,3) - slower)).
    kwargs :
        thr_cn: int, optional
            a threshold coordination number value. If provided, only atoms with coordination < thr_cn are considered for the GCN calculation
            (e.g. useful if you only want to consider surface atoms in your calculation)
        dbulk: float, optional
            Bulk distance for strained aGCN (default is None), has to be provided if strained is True

    Returns
    -------
    sites : ndarray
        sites (atomic coordinates. If thr_cn is enabled, only coordinates of atoms selected by the coordination condition)
    gcns : ndarray
        Values of the atop GCN.
    """
    neigh_list, coord_numbers = coordination_number(coords, cut_off, neigh_list=True, pbc=pbc, box=box)
    n_atoms = len(coord_numbers)
    agcn = []
    sites=[]

    thr_cn = kwargs.get('thr_cn', None)
    dbulk = kwargs.get('dbulk', None)

    if strained and dbulk is None:
        raise ValueError('Please provide bulk nn distance (dbulk) to compute the strained edition of this function!')

    for i, atom_neighbors in enumerate(neigh_list):
        if thr_cn is not None and coord_numbers[i] >= thr_cn:
            continue
        sites.append(coords[i])

        if strained:
            sgcn=0
            for nb in neigh_list[i]:
                for nnb in neigh_list[nb]:
                    if pbc:
                        d_nb_nnb = pbc_distance(coords[nb], coords[nnb], box)
                    else:
                        d_nb_nnb= np.linalg.norm(coords[nb] - coords[nnb])
                    sgcn += dbulk/d_nb_nnb
#            self_sgcn=0
#            for nb in neigh_list[i]:
#                break
#                d_nb_nnb= np.linalg.norm(coords[nb] - coords[i])
#                self_sgcn += dbulk/d_nb_nnb
#            agcn[i]=((sgcn-self_sgcn)/cn_max)
            agcn.append( sgcn/cn_max )
        else:
            agcn_i = sum(coord_numbers[neigh] for neigh in atom_neighbors)# - coord_numbers[i]
            agcn.append(agcn_i / cn_max)

    return np.array(sites), np.array(agcn)

bridge_gcn(coords, cut_off, thr_cn, dbulk=None, cn_max=18.0, phantom=True, strained=False, pbc=False, box=None)

Identifies bridge absorption sites and computes the Generalized Coordination Number (GCN) for a site. The GCN is defined as the sum of the coordination numbers of the neighbors of the two atoms forming the site, counted only once, divided by the reference maximum cn.

Parameters:

Name	Type	Description	Default
coords	ndarray	Array of the coordinates of the atoms forming the system.	required
cut_off	float	The cutoff distance for determining nearest neighbors.	required
cn_max	float	Maximum typical coordination number in the specific system (default is 18.0 - ok for FCC systems).	18.0
phantom	bool	If True, also returns the coordinates of the midpoints between pairs ('phantom' atoms indicating the bridge sites) (default is True).	True
thr_cn	int	a threshold coordination number value. Only atoms with coordination < thr_cn are considered for the GCN calculation (only surface bridge sites are considered)	required
strained	bool	if True, computes the strained aGCN.	False
dbulk	float	Bulk distance for strained aGCN (default is None), has to be provided if strained is True	None
pbc	bool	Whether to apply periodic boundary conditions. Defaults to False.	False
box	ndarray	Simulation box size (either [Lx, Ly, Lz] or [[xmin, xmax], [ymin, ymax], [zmin, zmax]] or 3 cell vectors (shape (3,3) - slower)).	None

Returns:

Name	Type	Description
sites	ndarray	Coordinates of the midpoints. Only returned if phantom is True
pairs	list	List of pairs.
bgcns	ndarray	Values of the bridge GCN ordered as the pairs.

Source code in snow/descriptors/coordination.py

def bridge_gcn(coords: np.ndarray, 
               cut_off: float, 
               thr_cn: int, 
               dbulk : float = None, 
               cn_max = 18.0,
               phantom: bool=True, 
               strained: bool = False, 
               pbc: bool = False, 
               box = None)-> tuple:
    """
    Identifies bridge absorption sites and computes the Generalized Coordination Number (GCN)
    for a site. The GCN is defined as the sum of the coordination numbers of the neighbors
    of the two atoms forming the site, counted only once, divided by the reference maximum cn.

    Parameters
    ----------
    coords : ndarray
        Array of the coordinates of the atoms forming the system.
    cut_off : float
        The cutoff distance for determining nearest neighbors.
    cn_max : float, default 18.0
        Maximum typical coordination number in the specific system (default is 18.0 - ok for FCC systems).
    phantom : bool, default True
        If True, also returns the coordinates of the midpoints between pairs ('phantom' atoms indicating the bridge sites)
        (default is True).
    thr_cn : int
        a threshold coordination number value. Only atoms with coordination < thr_cn are considered for the GCN calculation
        (only surface bridge sites are considered)
    strained : bool, default False
        if True, computes the strained aGCN.
    dbulk : float, optional
        Bulk distance for strained aGCN (default is None), has to be provided if strained is True
    pbc : bool, default False
        Whether to apply periodic boundary conditions. Defaults to False.
    box : np.ndarray, optional
        Simulation box size (either [Lx, Ly, Lz] or [[xmin, xmax], [ymin, ymax], [zmin, zmax]] or 3 cell vectors (shape (3,3) - slower)).

    Returns
    -------
    sites : ndarray
        Coordinates of the midpoints. Only returned if `phantom` is True
    pairs : list
        List of pairs.
    bgcns : ndarray
        Values of the bridge GCN ordered as the pairs.
    """

    #sanity check
    if strained and dbulk is None:
        raise ValueError('Please provide bulk nn distance (dbulk) to compute the strained edition of this function!')


    pairs = pair_list(coords=coords, cut_off=cut_off, pbc=pbc, box=box)
    neigh_list, coord_numb = coordination_number(coords=coords, cut_off=cut_off, neigh_list=True, pbc=pbc, box=box)
    #b_gcn = np.zeros(len(pairs))
    b_gcn=[]
    sites=[]
    for i, p in enumerate(pairs):
        if not (coord_numb[p[0]] < thr_cn and coord_numb[p[1]] < thr_cn):
            continue
        neigh_1 = neigh_list[p[0]]
        neigh_2 = neigh_list[p[1]]
        neigh_unique_12 = np.unique(np.concatenate((neigh_1, neigh_2)))
        if strained:
            sgcn = 0
            for nb in neigh_unique_12:
                for nnb in neigh_list[nb]:
                    if pbc:
                        d_nb_nnb = pbc_distance(coords[nb], coords[nnb], box)
                    else:
                        d_nb_nnb= np.linalg.norm(coords[nb] - coords[nnb])
                    sgcn += dbulk/d_nb_nnb
            self_sgcn=0
            for nb in p :
                for nnb in neigh_list[nb]:
                    if pbc:
                        d_nb_nnb = pbc_distance(coords[nb], coords[nnb], box)
                    else:
                        d_nb_nnb= np.linalg.norm(coords[nb] - coords[nnb])
                    self_sgcn += dbulk/d_nb_nnb
            b_gcn.append((sgcn-self_sgcn)/cn_max)
        else:
            b_gcn_i = sum(coord_numb[neigh] for neigh in neigh_unique_12) - (coord_numb[p[0]] + coord_numb[p[1]])
            b_gcn.append(b_gcn_i / cn_max)
        if phantom:
            pos_1 = coords[p[0]]
            pos_2 = coords[p[1]]
            sites.append((pos_1 + pos_2) / 2)
    if phantom:
        return np.asarray(sites), pairs, b_gcn
    else:
        return pairs, b_gcn

three_hollow_gcn(coords, cut_off, thr_cn, phantom=True, dbulk=None, cn_max=22.0, strained=False, pbc=None, box=None)

Finds the location of three-hollow sites and returns their location and GCN.

Parameters:

Name	Type	Description	Default
coords	ndarray	Array with the XYZ coordinates of the atoms, shape (n_atoms, 3).	required
cut_off	float	Cutoff distance for finding neighbors in angstrom.	required
thr_cn	int	a threshold coordination number value. Only atoms with coordination < thr_cn are considered for the GCN calculation (only surface hollow sites are considered).	required
phantom	bool	If True, also returns the coordinates of the midpoints between triplets ('phantom' atoms indicating the 3-hollow sites)	True
dbulk	float	Bulk distance for strained aGCN (default is None), has to be provided if strained is True	None
cn_max	float	Maximum typical coordination number in the specific system (default is 22.0 - ok for FCC materials).	22.0
strained	bool	if True, computes the strained aGCN.	False
pbc	bool	Whether to apply periodic boundary conditions.	False
box	ndarray	Simulation box size (either [Lx, Ly, Lz] or [[xmin, xmax], [ymin, ymax], [zmin, zmax]] or 3 cell vectors (shape (3,3) - slower)).	None

Returns:

Name	Type	Description
sites	ndarray	Coordinates of the triplets midpoints. Only returned if phantom is True
triplets	list[list]	List of atomic indexes labelling triplets.
th_gcns	ndarray	Values of the three-hollow GCN ordered as the pairs.

Source code in snow/descriptors/coordination.py

def three_hollow_gcn(coords: np.ndarray, 
                     cut_off: float, 
                     thr_cn: int, 
                     phantom: bool=True, 
                     dbulk: float = None, 
                     cn_max: float = 22.0, 
                     strained: bool = False, 
                     pbc: bool = None, 
                     box = None) -> tuple:
    """
    Finds the location of three-hollow sites and returns their location and GCN.

    Parameters
    ----------
    coords : np.ndarray
        Array with the XYZ coordinates of the atoms, shape (n_atoms, 3).
    cut_off : float
        Cutoff distance for finding neighbors in angstrom.
    thr_cn : int
        a threshold coordination number value. Only atoms with coordination < thr_cn are considered for the GCN calculation
        (only surface hollow sites are considered).
    phantom : bool, default True
        If True, also returns the coordinates of the midpoints between triplets ('phantom' atoms indicating the 3-hollow sites)
    dbulk : float, optional
        Bulk distance for strained aGCN (default is None), has to be provided if strained is True
    cn_max : float, default 22.0
        Maximum typical coordination number in the specific system (default is 22.0 - ok for FCC materials).
    strained : bool, default False
        if True, computes the strained aGCN.
    pbc : bool, default False
        Whether to apply periodic boundary conditions.
    box : np.ndarray, optional
        Simulation box size (either [Lx, Ly, Lz] or [[xmin, xmax], [ymin, ymax], [zmin, zmax]] or 3 cell vectors (shape (3,3) - slower)).

    Returns
    -------
    sites : ndarray
        Coordinates of the triplets midpoints. Only returned if `phantom` is True
    triplets : list[list]
        List of atomic indexes labelling triplets.
    th_gcns : ndarray
        Values of the three-hollow GCN ordered as the pairs.            

    """

    #sanity check
    if strained and dbulk is None:
        raise ValueError('Please provide bulk nn distance (dbulk) to compute the strained edition of this function!')

    triplets = []
    sites = []
    th_gcn = []
    pairs = pair_list(coords=coords, cut_off=cut_off, pbc=pbc, box=box)
    # neighbor list and coordination number not compatible!
    neigh_list, coord_numb = coordination_number(coords=coords, cut_off=cut_off, neigh_list=True, pbc=pbc, box=box)
    for i, p in enumerate(pairs):
        progress_bar(i, len(pairs) - 1, 50)
        # we check that both are at the surface
        # thr_cn is a threshold on cn to check if we are at surface
        if not (coord_numb[p[0]] < thr_cn and coord_numb[p[1]] < thr_cn):
            continue
        neigh_1 = neigh_list[p[0]]
        neigh_2 = neigh_list[p[1]]
        neigh_unique_12 = np.intersect1d(neigh_1, neigh_2, assume_unique=True)
        for cn in neigh_unique_12:
            if cn != p[0] and cn != p[1] and coord_numb[cn] < thr_cn:
                new_triplet = sorted([p[0], p[1], cn])
                if not new_triplet in triplets:  # this is a new surface triplet
                    triplets.append(new_triplet)
                    sites.append((coords[p[0]] + coords[p[1]] + coords[cn]) / 3)
                    neigh_unique_triplet = []
                    for idx in new_triplet:
                        neigh_unique_triplet += neigh_list[idx]
                    neigh_unique_triplet = np.unique(neigh_unique_triplet)
                    if strained:
                        sgcn = 0
                        for nb in neigh_unique_triplet:
                            for nnb in neigh_list[nb]:
                                if pbc:
                                    d_nb_nnb = pbc_distance(coords[nb], coords[nnb], box)
                                else:
                                    d_nb_nnb= np.linalg.norm(coords[nb] - coords[nnb])
                                sgcn += dbulk / d_nb_nnb
                        self_sgcn = 0
                        for nb in new_triplet:
                            for nnb in neigh_list[nb]:
                                if pbc:
                                    d_nb_nnb = pbc_distance(coords[nb], coords[nnb], box)
                                else:
                                    d_nb_nnb= np.linalg.norm(coords[nb] - coords[nnb])
                                self_sgcn += dbulk / d_nb_nnb
                        th_gcn.append((sgcn - self_sgcn) / cn_max)
                    else:
                        th_gcn_i = sum(coord_numb[neigh] for neigh in neigh_unique_triplet) - sum(
                            coord_numb[neigh] for neigh in new_triplet)
                        th_gcn.append(th_gcn_i / cn_max)

    if phantom:
        return np.asarray(sites), triplets, th_gcn
    else:
        return triplets, th_gcn

four_hollow_gcn(coords, cut_off, thr_cn, phantom=True, dbulk=None, cn_max=26.0, strained=False, pbc=False, box=None)

Finds the location of four-hollow sites and returns their location and GCN

Parameters:

Name	Type	Description	Default
coords	ndarray	Array with the XYZ coordinates of the atoms, shape (n_atoms, 3).	required
cut_off	float	Cutoff distance for finding neighbors in angstrom.	required
thr_cn	int	a threshold coordination number value. Only atoms with coordination < thr_cn are considered for the GCN calculation (e.g. only surface atoms are considered)	required
phantom	bool	If True, also returns the coordinates of the midpoints between fourplets ('phantom' atoms indicating the 4-hollow sites)	True
dbulk	float	Bulk distance for strained aGCN (default is None), has to be provided if strained is True	None
cn_max	float	Maximum typical coordination number in the specific system (default is 26.0 - ok for FCC materials).	26.0
strained	bool	if True, computes the strained aGCN.	False
pbc	bool	Whether to apply periodic boundary conditions. Defaults to False.	False
box	ndarray	Simulation box size (either [Lx, Ly, Lz] or [[xmin, xmax], [ymin, ymax], [zmin, zmax]] or 3 cell vectors (shape (3,3) - slower)).	None

Returns:

Name	Type	Description
	tuple
sites	ndarray	Coordinates of the fourplets midpoints. Only returned if phantom is True
fourplets	list[list]	List of atomic indexes labelling fourplets.
fh_gcns	ndarray	Values of the four-hollow GCN ordered as the pairs.

Source code in snow/descriptors/coordination.py

def four_hollow_gcn(coords: np.ndarray, 
                    cut_off: float, 
                    thr_cn: int, 
                    phantom: bool=True,
                    dbulk: float = None,
                    cn_max: float = 26.0, 
                    strained: bool = False, 
                    pbc: bool = False, 
                    box = None) -> tuple:
    """
    Finds the location of four-hollow sites and returns their location and GCN

    Parameters
    ----------
    coords: np.ndarray
        Array with the XYZ coordinates of the atoms, shape (n_atoms, 3).
    cut_off : float
        Cutoff distance for finding neighbors in angstrom.
    thr_cn : int
        a threshold coordination number value. Only atoms with coordination < thr_cn are considered for the GCN calculation
        (e.g. only surface atoms are considered)
    phantom : bool, default True
        If True, also returns the coordinates of the midpoints between fourplets ('phantom' atoms indicating the 4-hollow sites)
    dbulk: float, optional
        Bulk distance for strained aGCN (default is None), has to be provided if strained is True
    cn_max: float
        Maximum typical coordination number in the specific system (default is 26.0 - ok for FCC materials).
    strained : bool, default False
        if True, computes the strained aGCN.
    pbc : bool, default False
        Whether to apply periodic boundary conditions. Defaults to False.
    box : np.ndarray, optional
        Simulation box size (either [Lx, Ly, Lz] or [[xmin, xmax], [ymin, ymax], [zmin, zmax]] or 3 cell vectors (shape (3,3) - slower)).


    Returns
    -------
    tuple
    sites : ndarray
        Coordinates of the fourplets midpoints. Only returned if `phantom` is True
    fourplets : list[list]
        List of atomic indexes labelling fourplets.
    fh_gcns : ndarray
        Values of the four-hollow GCN ordered as the pairs.  

    """

    #sanity check
    if strained and dbulk is None:
        raise ValueError('Please provide bulk nn distance (dbulk) to compute the strained edition of this function!')


    fours = []
    sites = []
    fh_gcn = []
    pairs = pair_list(coords=coords, cut_off=cut_off, pbc=pbc, box=box)
    # neighbor list and coordination number not compatible!
    neigh_list, coord_numb = coordination_number(coords=coords, cut_off=cut_off, 
                                                 neigh_list=True, pbc=pbc, box=box)
    snb, _ = coordination_number(coords=coords, cut_off=cut_off * 1.3, neigh_list=True, pbc=pbc, box=box)
    current = 0
    for j in range(len(pairs)):
        for k in range(j):
            progress_bar(current, len(pairs) * (len(pairs) - 1) / 2)
            current += 1
            indices = [pairs[j][0], pairs[j][1], pairs[k][0], pairs[k][1]]
            check = False
            for idx in indices:  # check if all atoms are in surface
                if coord_numb[idx] >= thr_cn:
                    check = True
                # print(coord_numb[idx],thr_cn,check)
            if check:
                continue
            check = len(set(indices))
            if check != 4:
                continue  # got a shared atom
            # check if pairs are common atmostsecond neighbors
            common_nb_j = np.intersect1d(snb[indices[0]], snb[indices[1]], assume_unique=True)
            if (indices[2] in common_nb_j) and (indices[3] in common_nb_j):
                new_fours = sorted(indices)
                if not new_fours in fours:
                    fours.append(new_fours)
                    newsite = np.array([0., 0., 0.])
                    for at in new_fours:
                        newsite += coords[at]
                    sites.append(newsite / 4.0)
                    neigh_unique_four = []
                    for idx in new_fours:
                        neigh_unique_four += neigh_list[idx]
                    neigh_unique_four = np.unique(neigh_unique_four)
                    if strained:
                        sgcn = 0
                        for nb in neigh_unique_four:
                            for nnb in neigh_list[nb]:
                                if pbc:
                                    d_nb_nnb = pbc_distance(coords[nb], coords[nnb], box)
                                else:
                                    d_nb_nnb= np.linalg.norm(coords[nb] - coords[nnb])
                                sgcn += dbulk / d_nb_nnb
                        self_sgcn = 0
                        for nb in new_fours:
                            for nnb in neigh_list[nb]:
                                if pbc:
                                    d_nb_nnb = pbc_distance(coords[nb], coords[nnb], box)
                                else:
                                    d_nb_nnb= np.linalg.norm(coords[nb] - coords[nnb])
                                self_sgcn += dbulk / d_nb_nnb
                        fh_gcn.append((sgcn - self_sgcn) / cn_max)
                    else:
                        fh_gcn_i = sum(coord_numb[neigh] for neigh in neigh_unique_four) - sum(
                            coord_numb[neigh] for neigh in new_fours)
                        fh_gcn.append(fh_gcn_i / cn_max)

    #print("\nDone four hollow")
    if phantom:
        return np.asarray(sites), fours, fh_gcn
    else:
        return fours, fh_gcn