For each facet, plot the clean surface energy against the most
            stable binding energy.
        Args:
            delu_dict (Dict): Dictionary of the chemical potentials to be set as
                constant. Note the key should be a sympy Symbol object of the
                format: Symbol("delu_el") where el is the name of the element.
            delu_default (float): Default value for all unset chemical potentials
            plot_eads (bool): Option to plot the adsorption energy (binding
                energy multiplied by number of adsorbates) instead.
            annotate_monolayer (bool): Whether or not to label each data point
                with its monolayer (adsorbate density per unit primiitve area)
            JPERM2 (bool): Whether to plot surface energy in /m^2 (True) or
                eV/A^2 (False)

        Returns:
            (Plot): Plot of clean surface energy vs binding energy for
                all facets.

def BE_vs_clean_SE(self, delu_dict, delu_default=0, plot_eads=False,
                       annotate_monolayer=True, JPERM2=False):
        """
        For each facet, plot the clean surface energy against the most
            stable binding energy.
        Args:
            delu_dict (Dict): Dictionary of the chemical potentials to be set as
                constant. Note the key should be a sympy Symbol object of the
                format: Symbol("delu_el") where el is the name of the element.
            delu_default (float): Default value for all unset chemical potentials
            plot_eads (bool): Option to plot the adsorption energy (binding
                energy multiplied by number of adsorbates) instead.
            annotate_monolayer (bool): Whether or not to label each data point
                with its monolayer (adsorbate density per unit primiitve area)
            JPERM2 (bool): Whether to plot surface energy in /m^2 (True) or
                eV/A^2 (False)

        Returns:
            (Plot): Plot of clean surface energy vs binding energy for
                all facets.
        """

        plt = pretty_plot(width=8, height=7)
        for hkl in self.all_slab_entries.keys():
            for clean_entry in self.all_slab_entries[hkl].keys():
                all_delu_dict = self.set_all_variables(delu_dict, delu_default)
                if self.all_slab_entries[hkl][clean_entry]:
                    clean_se = self.as_coeffs_dict[clean_entry]
                    se = sub_chempots(clean_se, all_delu_dict)
                    for ads_entry in self.all_slab_entries[hkl][clean_entry]:
                        ml = ads_entry.get_monolayer
                        be = ads_entry.gibbs_binding_energy(eads=plot_eads)

                        # Now plot the surface energy vs binding energy
                        plt.scatter(se, be)
                        if annotate_monolayer:
                            plt.annotate("%.2f" % (ml), xy=[se, be],
                                         xytext=[se, be])

        plt.xlabel(r"Surface energy ($J/m^2$)") if JPERM2 \
            else plt.xlabel(r"Surface energy ($eV/\AA^2$)")
        plt.ylabel("Adsorption Energy (eV)") if plot_eads \
            else plt.ylabel("Binding Energy (eV)")
        plt.tight_layout()
        plt.xticks(rotation=60)

        return plt