Creates a dictionary where each entry is a key pointing to a
        chemical potential range where the surface of that entry is stable.
        Does so by enumerating through all possible solutions (intersect)
        for surface energies of a specific facet.

        Args:
            chempot_range ([max_chempot, min_chempot]): Range to consider the
                stability of the slabs.
            ref_delu (sympy Symbol): The range stability of each slab is based
                on the chempot range of this chempot. Should be a sympy Symbol
                object of the format: Symbol("delu_el") where el is the name of
                the element
            no_doped (bool): Consider stability of clean slabs only.
            no_clean (bool): Consider stability of doped slabs only.
            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.
            miller_index (list): Miller index for a specific facet to get a
                dictionary for.
            dmu_at_0 (bool): If True, if the surface energies corresponding to
                the chemical potential range is between a negative and positive
                value, the value is a list of three chemical potentials with the
                one in the center corresponding a surface energy of 0. Uselful
                in identifying unphysical ranges of surface energies and their
                chemical potential range.
            return_se_dict (bool): Whether or not to return the corresponding
                dictionary of surface energies

def stable_u_range_dict(self, chempot_range, ref_delu, no_doped=True,
                            no_clean=False, delu_dict={}, miller_index=(),
                            dmu_at_0=False, return_se_dict=False):
        """
        Creates a dictionary where each entry is a key pointing to a
        chemical potential range where the surface of that entry is stable.
        Does so by enumerating through all possible solutions (intersect)
        for surface energies of a specific facet.

        Args:
            chempot_range ([max_chempot, min_chempot]): Range to consider the
                stability of the slabs.
            ref_delu (sympy Symbol): The range stability of each slab is based
                on the chempot range of this chempot. Should be a sympy Symbol
                object of the format: Symbol("delu_el") where el is the name of
                the element
            no_doped (bool): Consider stability of clean slabs only.
            no_clean (bool): Consider stability of doped slabs only.
            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.
            miller_index (list): Miller index for a specific facet to get a
                dictionary for.
            dmu_at_0 (bool): If True, if the surface energies corresponding to
                the chemical potential range is between a negative and positive
                value, the value is a list of three chemical potentials with the
                one in the center corresponding a surface energy of 0. Uselful
                in identifying unphysical ranges of surface energies and their
                chemical potential range.
            return_se_dict (bool): Whether or not to return the corresponding
                dictionary of surface energies
        """

        chempot_range = sorted(chempot_range)
        stable_urange_dict, se_dict = {}, {}

        # Get all entries for a specific facet
        for hkl in self.all_slab_entries.keys():
            entries_in_hkl = []
            # Skip this facet if this is not the facet we want
            if miller_index and hkl != tuple(miller_index):
                continue
            if not no_clean:
                entries_in_hkl.extend([clean for clean in self.all_slab_entries[hkl]])
            if not no_doped:
                for entry in self.all_slab_entries[hkl]:
                    entries_in_hkl.extend([ads_entry for ads_entry in
                                           self.all_slab_entries[hkl][entry]])

            for entry in entries_in_hkl:
                stable_urange_dict[entry] = []
                se_dict[entry] = []
            # if there is only one entry for this facet, then just give it the
            # default urange, you can't make combinations with just 1 item
            if len(entries_in_hkl) == 1:
                stable_urange_dict[entries_in_hkl[0]] = chempot_range
                u1, u2 = delu_dict.copy(), delu_dict.copy()
                u1[ref_delu], u2[ref_delu] = chempot_range[0], chempot_range[1]
                se = self.as_coeffs_dict[entries_in_hkl[0]]
                se_dict[entries_in_hkl[0]] = [sub_chempots(se, u1), sub_chempots(se, u2)]
                continue

            for pair in itertools.combinations(entries_in_hkl, 2):
                # I'm assuming ref_delu was not set in delu_dict,
                # so the solution should be for ref_delu
                solution = self.get_surface_equilibrium(pair, delu_dict=delu_dict)

                # Check if this solution is stable
                if not solution:
                    continue
                new_delu_dict = delu_dict.copy()
                new_delu_dict[ref_delu] = solution[ref_delu]
                stable_entry, gamma = self.get_stable_entry_at_u(hkl, new_delu_dict,
                                                                 no_doped=no_doped,
                                                                 no_clean=no_clean)
                if stable_entry not in pair:
                    continue

                # Now check if the solution is within the chempot range
                if not (chempot_range[0] <= solution[ref_delu] <= chempot_range[1]):
                    continue

                for entry in pair:
                    stable_urange_dict[entry].append(solution[ref_delu])
                    se_dict[entry].append(gamma)

            # Now check if all entries have 2 chempot values. If only
            # one, we need to set the other value as either the upper
            # limit or lower limit of the user provided chempot_range
            new_delu_dict = delu_dict.copy()
            for u in chempot_range:
                new_delu_dict[ref_delu] = u
                entry, gamma = self.get_stable_entry_at_u(hkl, delu_dict=new_delu_dict,
                                                          no_doped=no_doped,
                                                          no_clean=no_clean)
                stable_urange_dict[entry].append(u)
                se_dict[entry].append(gamma)

        if dmu_at_0:
            for entry in se_dict.keys():
                # if se are of opposite sign, determine chempot when se=0.
                # Useful for finding a chempot range where se is unphysical
                if not stable_urange_dict[entry]:
                    continue
                if se_dict[entry][0] * se_dict[entry][1] < 0:
                    # solve for gamma=0
                    se = self.as_coeffs_dict[entry]
                    se_dict[entry].append(0)
                    stable_urange_dict[entry].append(solve(sub_chempots(se, delu_dict),
                                                           ref_delu)[0])

        # sort the chempot ranges for each facet
        for entry in stable_urange_dict.keys():
            se_dict[entry] = [se for i, se in sorted(zip(stable_urange_dict[entry],
                                                         se_dict[entry]))]
            stable_urange_dict[entry] = sorted(stable_urange_dict[entry])

        if return_se_dict:
            return stable_urange_dict, se_dict
        else:
            return stable_urange_dict