Function that helps visualize the slab in a 2-D plot, for
    convenient viewing of output of AdsorbateSiteFinder.

    Args:
        slab (slab): Slab object to be visualized
        ax (axes): matplotlib axes with which to visualize
        scale (float): radius scaling for sites
        repeat (int): number of repeating unit cells to visualize
        window (float): window for setting the axes limits, is essentially
            a fraction of the unit cell limits
        draw_unit_cell (bool): flag indicating whether or not to draw cell
        decay (float): how the alpha-value decays along the z-axis

def plot_slab(slab, ax, scale=0.8, repeat=5, window=1.5,
              draw_unit_cell=True, decay=0.2, adsorption_sites=True):
    """
    Function that helps visualize the slab in a 2-D plot, for
    convenient viewing of output of AdsorbateSiteFinder.

    Args:
        slab (slab): Slab object to be visualized
        ax (axes): matplotlib axes with which to visualize
        scale (float): radius scaling for sites
        repeat (int): number of repeating unit cells to visualize
        window (float): window for setting the axes limits, is essentially
            a fraction of the unit cell limits
        draw_unit_cell (bool): flag indicating whether or not to draw cell
        decay (float): how the alpha-value decays along the z-axis
    """
    orig_slab = slab.copy()
    slab = reorient_z(slab)
    orig_cell = slab.lattice.matrix.copy()
    if repeat:
        slab.make_supercell([repeat, repeat, 1])
    coords = np.array(sorted(slab.cart_coords, key=lambda x: x[2]))
    sites = sorted(slab.sites, key=lambda x: x.coords[2])
    alphas = 1 - decay * (np.max(coords[:, 2]) - coords[:, 2])
    alphas = alphas.clip(min=0)
    corner = [0, 0, slab.lattice.get_fractional_coords(coords[-1])[-1]]
    corner = slab.lattice.get_cartesian_coords(corner)[:2]
    verts = orig_cell[:2, :2]
    lattsum = verts[0] + verts[1]
    # Draw circles at sites and stack them accordingly
    for n, coord in enumerate(coords):
        r = sites[n].specie.atomic_radius * scale
        ax.add_patch(patches.Circle(coord[:2] - lattsum * (repeat // 2),
                                    r, color='w', zorder=2 * n))
        color = color_dict[sites[n].species_string]
        ax.add_patch(patches.Circle(coord[:2] - lattsum * (repeat // 2), r,
                                    facecolor=color, alpha=alphas[n],
                                    edgecolor='k', lw=0.3, zorder=2 * n + 1))
    # Adsorption sites
    if adsorption_sites:
        asf = AdsorbateSiteFinder(orig_slab)
        ads_sites = asf.find_adsorption_sites()['all']
        sop = get_rot(orig_slab)
        ads_sites = [sop.operate(ads_site)[:2].tolist()
                     for ads_site in ads_sites]
        ax.plot(*zip(*ads_sites), color='k', marker='x',
                markersize=10, mew=1, linestyle='', zorder=10000)
    # Draw unit cell
    if draw_unit_cell:
        verts = np.insert(verts, 1, lattsum, axis=0).tolist()
        verts += [[0., 0.]]
        verts = [[0., 0.]] + verts
        codes = [Path.MOVETO, Path.LINETO, Path.LINETO,
                 Path.LINETO, Path.CLOSEPOLY]
        verts = [(np.array(vert) + corner).tolist() for vert in verts]
        path = Path(verts, codes)
        patch = patches.PathPatch(path, facecolor='none', lw=2,
                                  alpha=0.5, zorder=2 * n + 2)
        ax.add_patch(patch)
    ax.set_aspect("equal")
    center = corner + lattsum / 2.
    extent = np.max(lattsum)
    lim_array = [center - extent * window, center + extent * window]
    x_lim = [ele[0] for ele in lim_array]
    y_lim = [ele[1] for ele in lim_array]
    ax.set_xlim(x_lim)
    ax.set_ylim(y_lim)
    return ax