def find_adsorption_sites(self, distance=2.0, put_inside=True,
                              symm_reduce=1e-2, near_reduce=1e-2,
                              positions=['ontop', 'bridge', 'hollow'],
                              no_obtuse_hollow=True):
        """
        Finds surface sites according to the above algorithm.  Returns
        a list of corresponding cartesian coordinates.

        Args:
            distance (float): distance from the coordinating ensemble
                of atoms along the miller index for the site (i. e.
                the distance from the slab itself)
            put_inside (bool): whether to put the site inside the cell
            symm_reduce (float): symm reduction threshold
            near_reduce (float): near reduction threshold
            positions (list): which positions to include in the site finding
                "ontop": sites on top of surface sites
                "bridge": sites at edges between surface sites in Delaunay
                    triangulation of surface sites in the miller plane
                "hollow": sites at centers of Delaunay triangulation faces
                "subsurface": subsurface positions projected into miller plane
            no_obtuse_hollow (bool): flag to indicate whether to include
                obtuse triangular ensembles in hollow sites
        """
        ads_sites = {k: [] for k in positions}
        if 'ontop' in positions:
            ads_sites['ontop'] = [s.coords for s in self.surface_sites]
        if 'subsurface' in positions:
            # Get highest site
            ref = self.slab.sites[np.argmax(self.slab.cart_coords[:, 2])]
            # Project diff between highest site and subs site into miller
            ss_sites = [self.mvec * np.dot(ref.coords - s.coords, self.mvec)
                        + s.coords for s in self.subsurface_sites()]
            ads_sites['subsurface'] = ss_sites
        if 'bridge' in positions or 'hollow' in positions:
            mesh = self.get_extended_surface_mesh()
            sop = get_rot(self.slab)
            dt = Delaunay([sop.operate(m.coords)[:2] for m in mesh])
            # TODO: refactor below to properly account for >3-fold
            for v in dt.simplices:
                if -1 not in v:
                    dots = []
                    for i_corner, i_opp in zip(range(3), ((1, 2), (0, 2), (0, 1))):
                        corner, opp = v[i_corner], [v[o] for o in i_opp]
                        vecs = [mesh[d].coords - mesh[corner].coords for d in opp]
                        vecs = [vec / np.linalg.norm(vec) for vec in vecs]
                        dots.append(np.dot(*vecs))
                        # Add bridge sites at midpoints of edges of D. Tri
                        if 'bridge' in positions:
                            ads_sites["bridge"].append(
                                self.ensemble_center(mesh, opp))
                    # Prevent addition of hollow sites in obtuse triangles
                    obtuse = no_obtuse_hollow and (np.array(dots) < 1e-5).any()
                    # Add hollow sites at centers of D. Tri faces
                    if 'hollow' in positions and not obtuse:
                        ads_sites['hollow'].append(
                            self.ensemble_center(mesh, v))
        ads_sites['all'] = sum(ads_sites.values(), [])
        for key, sites in ads_sites.items():
            # Pare off outer sites for bridge/hollow
            if key in ['bridge', 'hollow']:
                frac_coords = [self.slab.lattice.get_fractional_coords(ads_site)
                               for ads_site in sites]
                frac_coords = [frac_coord for frac_coord in frac_coords
                               if (frac_coord[0] > 1 and frac_coord[0] < 4
                                   and frac_coord[1] > 1 and frac_coord[1] < 4)]
                sites = [self.slab.lattice.get_cartesian_coords(frac_coord)
                         for frac_coord in frac_coords]
            if near_reduce:
                sites = self.near_reduce(sites, threshold=near_reduce)
            if put_inside:
                sites = [put_coord_inside(self.slab.lattice, coord)
                         for coord in sites]
            if symm_reduce:
                sites = self.symm_reduce(sites, threshold=symm_reduce)
            sites = [site + distance * self.mvec for site in sites]

            ads_sites[key] = sites
        return ads_sites