def monolayer_vs_BE(self, plot_eads=False):
        """
        Plots the binding energy energy as a function of monolayers (ML), i.e.
            the fractional area adsorbate density for all facets. For each
            facet at a specific monlayer, only plot the lowest binding energy.

        Args:
            plot_eads (bool): Option to plot the adsorption energy (binding
                 energy multiplied by number of adsorbates) instead.

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

        plt = pretty_plot(width=8, height=7)
        for hkl in self.all_slab_entries.keys():
            ml_be_dict = {}
            for clean_entry in self.all_slab_entries[hkl].keys():
                if self.all_slab_entries[hkl][clean_entry]:
                    for ads_entry in self.all_slab_entries[hkl][clean_entry]:
                        if ads_entry.get_monolayer not in ml_be_dict.keys():
                            ml_be_dict[ads_entry.get_monolayer] = 1000
                        be = ads_entry.gibbs_binding_energy(eads=plot_eads)
                        if be < ml_be_dict[ads_entry.get_monolayer]:
                            ml_be_dict[ads_entry.get_monolayer] = be
            # sort the binding energies and monolayers
            # in order to properly draw a line plot
            vals = sorted(ml_be_dict.items())
            monolayers, BEs = zip(*vals)
            plt.plot(monolayers, BEs, '-o',
                     c=self.color_dict[clean_entry], label=hkl)

        adsorbates = tuple(ads_entry.ads_entries_dict.keys())
        plt.xlabel(" %s" * len(adsorbates) % adsorbates + " Coverage (ML)")
        plt.ylabel("Adsorption Energy (eV)") if plot_eads \
            else plt.ylabel("Binding Energy (eV)")
        plt.legend()
        plt.tight_layout()

        return plt