def get_chempot_range_stability_phase(self, target_comp, open_elt):
        """
        returns a set of chemical potentials corresponding to the max and min
        chemical potential of the open element for a given composition. It is
        quite common to have for instance a ternary oxide (e.g., ABO3) for
        which you want to know what are the A and B chemical potential leading
        to the highest and lowest oxygen chemical potential (reducing and
        oxidizing conditions). This is useful for defect computations.

        Args:
            target_comp: A Composition object
            open_elt: Element that you want to constrain to be max or min

        Returns:
             {Element:(mu_min,mu_max)}: Chemical potentials are given in
             "absolute" values (i.e., not referenced to 0)
        """
        muref = np.array([self.el_refs[e].energy_per_atom
                          for e in self.elements if e != open_elt])
        chempot_ranges = self.get_chempot_range_map(
            [e for e in self.elements if e != open_elt])
        for e in self.elements:
            if not e in target_comp.elements:
                target_comp = target_comp + Composition({e: 0.0})
        coeff = [-target_comp[e] for e in self.elements if e != open_elt]
        max_open = -float('inf')
        min_open = float('inf')
        max_mus = None
        min_mus = None
        for e in chempot_ranges.keys():
            if e.composition.reduced_composition == \
                    target_comp.reduced_composition:
                multiplicator = e.composition[open_elt] / target_comp[open_elt]
                ef = e.energy / multiplicator
                all_coords = []
                for s in chempot_ranges[e]:
                    for v in s._coords:
                        all_coords.append(v)
                        if (np.dot(v + muref, coeff) + ef) / target_comp[
                                open_elt] > max_open:
                            max_open = (np.dot(v + muref, coeff) + ef) / \
                                       target_comp[open_elt]
                            max_mus = v
                        if (np.dot(v + muref, coeff) + ef) / target_comp[
                                open_elt] < min_open:
                            min_open = (np.dot(v + muref, coeff) + ef) / \
                                       target_comp[open_elt]
                            min_mus = v
        elts = [e for e in self.elements if e != open_elt]
        res = {}
        for i in range(len(elts)):
            res[elts[i]] = (min_mus[i] + muref[i], max_mus[i] + muref[i])
        res[open_elt] = (min_open, max_open)
        return res