Finds bonded atoms and returns a adjacency matrix of bonded atoms.

    Author: "Gowoon Cheon"
    Email: "gcheon@stanford.edu"

    Args:
        struct (Structure): Input structure
        tolerance: length in angstroms used in finding bonded atoms. Two atoms
            are considered bonded if (radius of atom 1) + (radius of atom 2) +
            (tolerance) < (distance between atoms 1 and 2). Default
            value = 0.45, the value used by JMol and Cheon et al.
        ldict: dictionary of bond lengths used in finding bonded atoms. Values
            from JMol are used as default

    Returns:
        (np.ndarray): A numpy array of shape (number of atoms, number of atoms);
        If any image of atom j is bonded to atom i with periodic boundary
        conditions, the matrix element [atom i, atom j] is 1.

def find_connected_atoms(struct, tolerance=0.45, ldict=JmolNN().el_radius):
    """
    Finds bonded atoms and returns a adjacency matrix of bonded atoms.

    Author: "Gowoon Cheon"
    Email: "gcheon@stanford.edu"

    Args:
        struct (Structure): Input structure
        tolerance: length in angstroms used in finding bonded atoms. Two atoms
            are considered bonded if (radius of atom 1) + (radius of atom 2) +
            (tolerance) < (distance between atoms 1 and 2). Default
            value = 0.45, the value used by JMol and Cheon et al.
        ldict: dictionary of bond lengths used in finding bonded atoms. Values
            from JMol are used as default

    Returns:
        (np.ndarray): A numpy array of shape (number of atoms, number of atoms);
        If any image of atom j is bonded to atom i with periodic boundary
        conditions, the matrix element [atom i, atom j] is 1.
    """
    n_atoms = len(struct.species)
    fc = np.array(struct.frac_coords)
    fc_copy = np.repeat(fc[:, :, np.newaxis], 27, axis=2)
    neighbors = np.array(list(itertools.product([0, 1, -1], [0, 1, -1], [0, 1, -1]))).T
    neighbors = np.repeat(neighbors[np.newaxis, :, :], 1, axis=0)
    fc_diff = fc_copy - neighbors
    species = list(map(str, struct.species))
    # in case of charged species
    for i, item in enumerate(species):
        if not item in ldict.keys():
            species[i] = str(Specie.from_string(item).element)
    latmat = struct.lattice.matrix
    connected_matrix = np.zeros((n_atoms,n_atoms))

    for i in range(n_atoms):
        for j in range(i + 1, n_atoms):
            max_bond_length = ldict[species[i]] + ldict[species[j]] + tolerance
            frac_diff = fc_diff[j] - fc_copy[i]
            distance_ij = np.dot(latmat.T, frac_diff)
            # print(np.linalg.norm(distance_ij,axis=0))
            if sum(np.linalg.norm(distance_ij, axis=0) < max_bond_length) > 0:
                connected_matrix[i, j] = 1
                connected_matrix[j, i] = 1
    return connected_matrix