Returns the surface sites and their indices in a dictionary. The
        oriented unit cell of the slab will determine the coordination number
        of a typical site. We use VoronoiNN to determine the
        coordination number of bulk sites and slab sites. Due to the
        pathological error resulting from some surface sites in the
        VoronoiNN, we assume any site that has this error is a surface
        site as well. This will work for elemental systems only for now. Useful
        for analysis involving broken bonds and for finding adsorption sites.

            Args:
                tag (bool): Option to adds site attribute "is_surfsite" (bool)
                    to all sites of slab. Defaults to False

            Returns:
                A dictionary grouping sites on top and bottom of the slab
                together.
                {"top": [sites with indices], "bottom": [sites with indices}

        TODO:
            Is there a way to determine site equivalence between sites in a slab
            and bulk system? This would allow us get the coordination number of
            a specific site for multi-elemental systems or systems with more
            than one unequivalent site. This will allow us to use this for
            compound systems.

def get_surface_sites(self, tag=False):
        """
        Returns the surface sites and their indices in a dictionary. The
        oriented unit cell of the slab will determine the coordination number
        of a typical site. We use VoronoiNN to determine the
        coordination number of bulk sites and slab sites. Due to the
        pathological error resulting from some surface sites in the
        VoronoiNN, we assume any site that has this error is a surface
        site as well. This will work for elemental systems only for now. Useful
        for analysis involving broken bonds and for finding adsorption sites.

            Args:
                tag (bool): Option to adds site attribute "is_surfsite" (bool)
                    to all sites of slab. Defaults to False

            Returns:
                A dictionary grouping sites on top and bottom of the slab
                together.
                {"top": [sites with indices], "bottom": [sites with indices}

        TODO:
            Is there a way to determine site equivalence between sites in a slab
            and bulk system? This would allow us get the coordination number of
            a specific site for multi-elemental systems or systems with more
            than one unequivalent site. This will allow us to use this for
            compound systems.
        """

        from pymatgen.analysis.local_env import VoronoiNN

        # Get a dictionary of coordination numbers
        # for each distinct site in the structure
        a = SpacegroupAnalyzer(self.oriented_unit_cell)
        ucell = a.get_symmetrized_structure()
        cn_dict = {}
        v = VoronoiNN()
        unique_indices = [equ[0] for equ in ucell.equivalent_indices]

        for i in unique_indices:
            el = ucell[i].species_string
            if el not in cn_dict.keys():
                cn_dict[el] = []
            # Since this will get the cn as a result of the weighted polyhedra, the
            # slightest difference in cn will indicate a different environment for a
            # species, eg. bond distance of each neighbor or neighbor species. The
            # decimal place to get some cn to be equal.
            cn = v.get_cn(ucell, i, use_weights=True)
            cn = float('%.5f' %(round(cn, 5)))
            if cn not in cn_dict[el]:
                cn_dict[el].append(cn)

        v = VoronoiNN()

        surf_sites_dict, properties = {"top": [], "bottom": []}, []
        for i, site in enumerate(self):
            # Determine if site is closer to the top or bottom of the slab
            top = True if site.frac_coords[2] > self.center_of_mass[2] else False

            try:
                # A site is a surface site, if its environment does
                # not fit the environment of other sites
                cn = float('%.5f' %(round(v.get_cn(self, i, use_weights=True), 5)))
                if cn < min(cn_dict[site.species_string]):
                    properties.append(True)
                    key = "top" if top else "bottom"
                    surf_sites_dict[key].append([site, i])
                else:
                    properties.append(False)
            except RuntimeError:
                # or if pathological error is returned, indicating a surface site
                properties.append(True)
                key = "top" if top else "bottom"
                surf_sites_dict[key].append([site, i])

        if tag:
            self.add_site_property("is_surf_site", properties)
        return surf_sites_dict