def get_site_dict(self, ads_pos):
        """Get dictionary with high symmetry sites close to adsorbate
        position. 
        Top sites: Optained from the atomic positions of
                   the top layer of the slab. 
        Bridge sites: The average position of atomic pairs
        Hollow sites: Optained as the Voronoi vertices
        4-fold sites: Assigned when Voronoi vertives overlap
        """

        # Get top layer
        C = self.B[-self.ntop - 1:-1]
        SC = C * (3, 3, 1)

        D = {}
        cell = C.get_cell()
        top_dist = 1

        """Top sites: atomic positions """
        for i, atom in \
            enumerate([atom for atom in SC if
                       np.linalg.norm(atom.position[:2] - ads_pos)
                       < top_dist]):
            k = 'top_site-{}'.format(i)
            D[k] = {}
            D[k]['pos'] = atom.position
            D[k]['sym'] = atom.symbol

        """Bridge sites: bewteen atomic pairs """
        bridge_dist = 0.5 * cell[0][0]
        b_index = 0
        for atomi in [atom for atom in SC if
                      np.linalg.norm(atom.position[:2] - ads_pos)
                      < bridge_dist]:
            pos1 = atomi.position
            # Atom bairs should be close to each other and close to
            # adsorbate
            for atom in \
                [atom for atom in SC if np.linalg.norm(atom.position - pos1)
                 < 1.5 * bridge_dist and
                 np.linalg.norm(atom.position[:2] - ads_pos) < bridge_dist
                 and not (atom.position == pos1).all()]:
                pos2 = atom.position
                # bridge site is average position
                bridge_pos = 0.5 * (pos1 + pos2)
                k = 'bridge_site-{}'.format(b_index)
                D[k] = {}
                D[k]['pos'] = bridge_pos
                D[k]['sym'] = atomi.symbol + '_' + atom.symbol
                b_index += 1

        """Hollow sites: Voronoi vertices """
        hollow_dist = 1
        vor = Voronoi(SC.positions[:, :2])
        vertices = vor.vertices
        # map vertices close to adsorbate
        close_v = [v for v in vertices if np.linalg.norm(v - ads_pos[:2])
                   < hollow_dist]
        h_index = 1
        ff_index = 1
        for v in close_v:
            # Check if vertices overlap
            close_v_v = [v0 for v0 in close_v if np.linalg.norm(v0 - v) < 0.5]
            if len(close_v_v) > 1:
                v_mean = np.mean(close_v_v, axis=0)
                k = '4fold_{}'.format(ff_index)
                D[k] = {}
                D[k]['pos'] = v_mean

                # Delete bridge sites overlapping with 4fold
                for key in [key for key in list(D.keys()) if 'bridge' in key]:
                    bridge_pos = D[key]['pos']
                    if np.linalg.norm(bridge_pos[:2] - v) < 0.3:
                        del D[key]

                ffold_max_dist = sorted([np.linalg.norm(v_mean - m[:2])
                                         for m in SC.positions])[4]
                symb = [atom.symbol for atom in SC if
                        np.linalg.norm(v_mean - atom.position[:2])
                        < ffold_max_dist]

                D[k]['sym'] = '_'.join(symb)
                ff_index += 1
            else:  # Regular hollow site
                k = 'hollow_site-{}'.format(h_index)
                D[k] = {}
                D[k]['pos'] = v
                hollow_max_dist = sorted([np.linalg.norm(v - m[:2])
                                          for m in SC.positions])[3]
                symb = [atom.symbol for atom in SC
                        if np.linalg.norm(v - atom.position[:2])
                        < hollow_max_dist]
                D[k]['sym'] = '_'.join(symb)
                h_index += 1

        return D