code: def monolayer_vs_BE(self, plot_eads=False):
          
  
          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
docs: 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.