<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 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