Adapted from the get_chempot_range_map() method in the PhaseDiagram
            class. Plot the chemical potential range map based on surface
            energy stability. Currently works only for 2-component PDs. At
            the moment uses a brute force method by enumerating through the
            range of the first element chempot with a specified increment
            and determines the chempot rangeo fht e second element for each
            SlabEntry. Future implementation will determine the chempot range
            map first by solving systems of equations up to 3 instead of 2.
        Args:
            elements (list): Sequence of elements to be considered as independent
                variables. E.g., if you want to show the stability ranges of
                all Li-Co-O phases wrt to duLi and duO, you will supply
                [Element("Li"), Element("O")]
            miller_index ([h, k, l]): Miller index of the surface we are interested in
            ranges ([[range1], [range2]]): List of chempot ranges (max and min values)
                for the first and second element.
            incr (int): Number of points to sample along the range of the first chempot
            no_doped (bool): Whether or not to include doped systems.
            no_clean (bool): Whether or not to include clean systems.
            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.
            annotate (bool): Whether to annotate each "phase" with the label of
                the entry. If no label, uses the reduced formula
            show_unphyiscal_only (bool): Whether to only show the shaded region where
                surface energy is negative. Useful for drawing other chempot range maps.

def surface_chempot_range_map(self, elements, miller_index, ranges,
                                  incr=50, no_doped=False, no_clean=False,
                                  delu_dict=None, plt=None, annotate=True,
                                  show_unphyiscal_only=False, fontsize=10):
        """
        Adapted from the get_chempot_range_map() method in the PhaseDiagram
            class. Plot the chemical potential range map based on surface
            energy stability. Currently works only for 2-component PDs. At
            the moment uses a brute force method by enumerating through the
            range of the first element chempot with a specified increment
            and determines the chempot rangeo fht e second element for each
            SlabEntry. Future implementation will determine the chempot range
            map first by solving systems of equations up to 3 instead of 2.
        Args:
            elements (list): Sequence of elements to be considered as independent
                variables. E.g., if you want to show the stability ranges of
                all Li-Co-O phases wrt to duLi and duO, you will supply
                [Element("Li"), Element("O")]
            miller_index ([h, k, l]): Miller index of the surface we are interested in
            ranges ([[range1], [range2]]): List of chempot ranges (max and min values)
                for the first and second element.
            incr (int): Number of points to sample along the range of the first chempot
            no_doped (bool): Whether or not to include doped systems.
            no_clean (bool): Whether or not to include clean systems.
            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.
            annotate (bool): Whether to annotate each "phase" with the label of
                the entry. If no label, uses the reduced formula
            show_unphyiscal_only (bool): Whether to only show the shaded region where
                surface energy is negative. Useful for drawing other chempot range maps.
        """

        # Set up
        delu_dict = delu_dict if delu_dict else {}
        plt = pretty_plot(12, 8) if not plt else plt
        el1, el2 = str(elements[0]), str(elements[1])
        delu1 = Symbol("delu_%s" % (str(elements[0])))
        delu2 = Symbol("delu_%s" % (str(elements[1])))
        range1 = ranges[0]
        range2 = ranges[1]

        # Find a range map for each entry (surface). This part is very slow, will
        # need to implement a more sophisticated method of getting the range map
        vertices_dict = {}
        for dmu1 in np.linspace(range1[0], range1[1], incr):
            # Get chemical potential range of dmu2 for each increment of dmu1
            new_delu_dict = delu_dict.copy()
            new_delu_dict[delu1] = dmu1
            range_dict, se_dict = self.stable_u_range_dict(range2, delu2, dmu_at_0=True,
                                                           miller_index=miller_index,
                                                           no_doped=no_doped,
                                                           no_clean=no_clean,
                                                           delu_dict=new_delu_dict,
                                                           return_se_dict=True)

            # Save the chempot range for dmu1 and dmu2
            for entry in range_dict.keys():
                if not range_dict[entry]:
                    continue
                if entry not in vertices_dict.keys():
                    vertices_dict[entry] = []

                selist = se_dict[entry]
                vertices_dict[entry].append({delu1: dmu1, delu2: [range_dict[entry], selist]})

        # Plot the edges of the phases
        for entry in vertices_dict.keys():

            xvals, yvals = [], []

            # Plot each edge of a phase within the borders
            for ii, pt1 in enumerate(vertices_dict[entry]):

                # Determine if the surface energy at this lower range
                # of dmu2 is negative. If so, shade this region.
                if len(pt1[delu2][1]) == 3:
                    if pt1[delu2][1][0] < 0:
                        neg_dmu_range = [pt1[delu2][0][0], pt1[delu2][0][1]]
                    else:
                        neg_dmu_range = [pt1[delu2][0][1], pt1[delu2][0][2]]
                    # Shade the threshold and region at which se<=0
                    plt.plot([pt1[delu1], pt1[delu1]], neg_dmu_range, 'k--')
                elif pt1[delu2][1][0] < 0 and pt1[delu2][1][1] < 0:
                    # Any chempot at at this point will result
                    # in se<0, shade the entire y range
                    if not show_unphyiscal_only:
                        plt.plot([pt1[delu1], pt1[delu1]], range2, 'k--')

                if ii == len(vertices_dict[entry]) - 1:
                    break
                pt2 = vertices_dict[entry][ii + 1]
                if not show_unphyiscal_only:
                    plt.plot([pt1[delu1], pt2[delu1]], [pt1[delu2][0][0], pt2[delu2][0][0]], 'k')

                # Need these values to get a good position for labelling phases
                xvals.extend([pt1[delu1], pt2[delu1]])
                yvals.extend([pt1[delu2][0][0], pt2[delu2][0][0]])

            # Plot the edge along the max x value
            pt = vertices_dict[entry][-1]
            delu1, delu2 = pt.keys()
            xvals.extend([pt[delu1], pt[delu1]])
            yvals.extend(pt[delu2][0])
            if not show_unphyiscal_only:
                plt.plot([pt[delu1], pt[delu1]], [pt[delu2][0][0], pt[delu2][0][-1]], 'k')

            if annotate:
                # Label the phases
                x = np.mean([max(xvals), min(xvals)])
                y = np.mean([max(yvals), min(yvals)])
                label = entry.label if entry.label else entry.composition.reduced_formula
                plt.annotate(label, xy=[x, y], xytext=[x, y], fontsize=fontsize)

        # Label plot
        plt.xlim(range1)
        plt.ylim(range2)
        plt.xlabel(r"$\Delta\mu_{%s} (eV)$" % (el1), fontsize=25)
        plt.ylabel(r"$\Delta\mu_{%s} (eV)$" % (el2), fontsize=25)
        plt.xticks(rotation=60)

        return plt