code: def cluster_sites(mol, tol, give_only_index=False):
      
      # Cluster works for dim > 2 data. We just add a dummy 0 for second
      # coordinate.
      dists = [[np.linalg.norm(site.coords), 0] for site in mol]
      import scipy.cluster as spcluster
      f = spcluster.hierarchy.fclusterdata(dists, tol, criterion='distance')
      clustered_dists = defaultdict(list)
      for i, site in enumerate(mol):
          clustered_dists[f[i]].append(dists[i])
      avg_dist = {label: np.mean(val) for label, val in clustered_dists.items()}
      clustered_sites = defaultdict(list)
      origin_site = None
      for i, site in enumerate(mol):
          if avg_dist[f[i]] < tol:
              if give_only_index:
                  origin_site = i
              else:
                  origin_site = site
          else:
              if give_only_index:
                  clustered_sites[
                      (avg_dist[f[i]], site.species)].append(i)
              else:
                  clustered_sites[
                      (avg_dist[f[i]], site.species)].append(site)
      return origin_site, clustered_sites
docs: Cluster sites based on distance and species type.
  
      Args:
          mol (Molecule): Molecule **with origin at center of mass**.
          tol (float): Tolerance to use.
  
      Returns:
          (origin_site, clustered_sites): origin_site is a site at the center
          of mass (None if there are no origin atoms). clustered_sites is a
          dict of {(avg_dist, species_and_occu): [list of sites]}